1 ## Copyright (C) 1995-2012 Kurt Hornik
3 ## This file is part of Octave.
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 3 of the License, or (at
8 ## your option) any later version.
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
20 ## @deftypefn {Function File} {} duplication_matrix (@var{n})
21 ## Return the duplication matrix
28 ## which is the unique
30 ## $n^2 \times n(n+1)/2$
33 ## @math{n^2} by @math{n*(n+1)/2}
37 ## $D_n * {\rm vech} (A) = {\rm vec} (A)$
40 ## @math{Dn vech (A) = vec (A)}
47 ## @math{n} by @math{n}
57 ## See Magnus and Neudecker (1988), Matrix differential calculus with
58 ## applications in statistics and econometrics.
61 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
62 ## Created: 8 May 1995
65 function d = duplication_matrix (n)
71 if (! (isscalar (n) && n > 0 && n == fix (n)))
72 error ("duplication_matrix: N must be a positive integer");
75 d = zeros (n * n, n * (n + 1) / 2);
77 ## It is clearly possible to make this a LOT faster!
80 d ((j - 1) * n + j, count + j) = 1;
82 d ((j - 1) * n + i, count + i) = 1;
83 d ((i - 1) * n + j, count + i) = 1;
85 count = count + n - j;
95 %! D = duplication_matrix (N);
96 %! assert (D * vech (B), vec (B), 1e-6);
97 %! assert (D * vech (C), vec (C), 1e-6);
104 %! D = duplication_matrix (N);
105 %! assert (D * vech (B), vec (B), 1e-6);
106 %! assert (D * vech (C), vec (C), 1e-6);
113 %! D = duplication_matrix (N);
114 %! assert (D * vech (B), vec (B), 1e-6);
115 %! assert (D * vech (C), vec (C), 1e-6);
117 %!error duplication_matrix ();
118 %!error duplication_matrix (0.5);
119 %!error duplication_matrix (-1);
120 %!error duplication_matrix (ones(1,4));