vech 2nd argument

[Rolf Fabian]

A second argument to 'vech.m' is of big value
at least for my personal linear algebra programs.

I added it, changed the help part nearly completely
and added a bunch of input checks. The previous
version had essentially none of them except
checking that matrix is square.

I don't see any reason not to implement a
second argument as 'vech(x,k)' into vech
also in order to emphasize its affinity to
'tril(x,k)' more clearly.

Rolf Fabian
< r fabian at jacobs-university dot de >



An example session:
octave-2.9.19.exe:15> x=reshape(1:16,4,4)
x =
    1    5    9   13
    2    6   10   14
    3    7   11   15
    4    8   12   16

octave-2.9.19.exe:16> vech(x,-1)
ans =
    2
    3
    4
    7
    8
   12

octave-2.9.19.exe:17> vech(x,-2)
ans =
   3
   4
   8

octave-2.9.19.exe:18> vech(x,-3)
ans =  4
octave-2.9.19.exe:19> vech(x,-4)
ans = [](0x1)
octave-2.9.19.exe:20> vech(x,-5)
error: vech: abs(k) cannot exceed dimension of square matrix.
error: ....
octave-2.9.19.exe:20> vech(x,1)
error: vech: case k > 0 not implemented.
error: ....


-------------------- snip -------------------
## Copyright (C) 1995, 1996, 1997, 1999, 2000, 2002, 2005, 2006, 2007
##               Kurt Hornik , (C) Rolf Fabian 2007
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {} vech (@var{x},@var{k})
## Return the column vector obtained by considering only lower triangular
## part of order @var{k} ( i.e. @code{tril(@var{x},@var{k})} ) elements
## of square matrix @var{x} and stacking them each column above the other.
## A single argument @var{k} implies @var{k}=0. Note that any nonzero
## @var{k} needs to be a negative integer with magnitude up to matrix
## dimension of @var{x}.
## @end deftypefn

## See Magnus and Neudecker (1988), Matrix differential calculus with
## applications in statistics and econometrics.

## Author KH <[e-mail]>
## Created: 8 May 1995
## Adapted-By: jwe
##
## Rolf Fabian < r dot fabian at jacobs-university dot de>
## Last modification : 2007-12-18 (based on V2.9.19)

function v = vech (x,k)

  if  (nargin < 1 || nargin > 2)
      print_usage ();
  endif
 
  if  ( nargin == 1 )
      k = 0;
  endif
 
  if  (  ndims (x) > 2 || isstruct (x) || iscell (x)
      || ndims (k) > 2 || isstruct (k) || iscell (k)
      || iscomplex (k) || ! isscalar (k) || k - round(k) )
      error ("vech: type mismatch of input.");
  endif

  if  (! (n = issquare (x)) )
      error ("vech: x must be square");
  endif

  if ( k>0 )
      error ("vech: case k > 0 not implemented.");
  elseif ( -k > n )
      error ("vech: abs(k) cannot exceed dimension of square matrix.");
  endif
   
  x( 1:(-k), : ) = [];
  ## delete uppermost abs(k) rows
  x( :, n+k+1:n ) = [];
  ## delete most right abs(k) columns

  ## This should be quicker than having an inner `for' loop as well.
  ## Ideally, vech should be written in C++.
  n = rows (x);
  v = zeros ((n+1)*n/2, 1);
  count = 0;
  for j = 1 : n
      i = j : n;
      v (count + i) = x (i, j);
      count = count + n - j;
  endfor

endfunction