LinAlg: matrix-normality check proposed, g-circulant matrices

[Rolf Fabian]

Most of you probably know that hermitian matrices have an
orthogonal system of eigenvectors. As a consequence
they may be always diagonalized via an unitary similarity
transformation.

>From a numerical point of view checking for matrix normality
prior to a matrix similarity transformation helps to avoid the
time-consuming calculation of an inverse matrix via 'inv',
which is needed to perform a 'general' similarity transformation.

Of course, a check for square matrix normality is required because
Octave hasn't one implemented. Thus I propose to implement my
'isnormal.m' function below.

Many of you most probably don't know that 'hermitian' matrices
correspond only to a very limited subset of 'normal' matrices.
Description file 'normal.txt' below lists about * 30 more classes * 
of 'normal' matrices which are all recognized by 'isnormal.m'.

g-circulant matrices are probably not that well-known.
I've added 'circular.m' in order to allow you to create them
for checking 'isnormal'. In my eyes, 'circular.m' is also worth
to be implemented into Octave (or Octave-forge).
Its 3rd argument is needed because some authors in the
corresponding literature define circulation along columns,
whereas others prefer circulation along rows.

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


_________________________________________________________

## Copyright (C) 1997 - 2008  Rolf Fabian
##
## This file is part of Octave.
##
## It is published under the 'GNU General Public License' of the 'Free
## Software Foundation', either Version 3 of the license, or any later
## version. This work is free software. You can redistribute it and/or
## modify it under the terms of the specified license.
##
## This work 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.
##
## In case you have not received a copy of the specified license along
## with this file, consult GNU-website <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} address@hidden,@var{m},@var{nfro}]} = isnormal
(@var{x}, @var{tol}).
## Return its dimension @var{n} if square matrix @var{x} is found
## to be @code{normal} within @var{tol}, otherwise 0. By default,
## @address@hidden sqrt(2) * @var{n} * eps} is used. Tolerance feature
## may be switched off completely by using @address@hidden
##
## Square matrix @var{x} pocesses @code{normal} property if it commutes
## with its own conjugate complex transposed matrix @var{x}'. Then its
## eigenvectors span complete address@hidden@var{n}] vector space and @var{x}
## is always @code{unitarily diagonalizable}.
##
## Actually, the normality check is performed by evaluation of
## @code{m = norm (@var{y} * @var{y}' - @var{y}' * @var{y}, 'fro')} where
## @address@hidden = @var{x} / @var{nfro}} corresponds to the input
## matrix scaled by its Frobenius norm. Scaling avoids over/underflow
## for very large/small input matrix norms during computation of
## implied matrix products.
## @seealso{circulant}
## @end deftypefn

## Author : Rolf Fabian <r dot fabian at Jacobs-University dot de>
## Last revision : 2008-02-19


function [ n, m, nfro ] = isnormal (x, tol)

if (nargin <1 || nargin >2)
   usage ("[n, m, nfo] = isnormal (x, tol)");
endif

m = 0;
n = issquare (x);

if (! n)
   error ("isnormal: square matrix input required.");
elseif (n == 1)
   ## scalars are always 'normal'
   nfro = abs (x);
   return;
endif

if (nargin < 2)
   ## default tolerance
   tol = sqrt (2) .* n .* eps;
endif

if (iscomplex (tol) || ! isscalar (tol) || tol <0 )
   error ("isnormal: argument 'tol' needs to be a positive (real) scalar.");
endif

nfro = norm (x,'fro');
## get Frobenius norm to be used as scale factor

if (nfro && (fx = finite (nfro)))
   x ./= nfro;
elseif (! nfro)
   ## detected x=zeros(n) matrix, which is always normal
   return;
endif

m = norm ( x * x' .- x' * x ,'fro');

if (m > tol || ! fx )
   ## force also n=0 (false) if x is not finite
   n = 0;
endif

endfunction

_________________________________________________________

## Author : Rolf Fabian <r dot fabian at jacobs-university dot de>
## Latest revision : 2008-02-19

The table below shows the occurence of the 'normal' property for
classes of square [nxn] matrices pocessing special properties
(like e.g. special symmetry patterns) in dependence of the
matrix complexity status (real, complex, imaginary).

The eigenvectors of 'normal' matrices always form an orthogonal
basis spanning complete [nxn] vector space. Thus, such
matrices are always unitarily diagonalizable.


# 1 : square mat is generally 'normal' (orthogonal eigensystem)
# 0 : square mat is generally  NOT 'normal' 

symm.
flag       real    complex    imaginary               matrix attribute
.........................................................................................
## standard symmetry patterns
's'          1           0            1                        'symmetric'
'sh'        1            1           1                        'hermitian'
'sh-'       1           1            1                       
'skewhermitian'
's-'         1           0            1                       
'skewsymmetric'

## Toeplitz symmetry patterns
't'          1            0           1                       
'toeplitzsymmetric'
'th'        1            1           1                       
'toeplitzhermitian'
'th-'       1            1           1                       
'toeplitzskewhermitian'
't-'         1            0           1                       
'toeplitzskewsymmetric'

## Hankel symmetry patterns
'h'          1           0           1                       
'hankelsymmetric'
'hh'        1           0            1                       
'hankelhermitian'
'hh-'       1           0            1                       
'hankelskewhermitian'
'h-'         1           0           1                       
'hankelskewsymmetric'

## Centrosymmetry class
'c*'         0          0          0                          all 4 patterns
of 'centrosymmetry' class

## Orthogonal                                            "x.' * x =
eye(dim)"
              1          0          1

## Unitary                                                 "x'  * x =
eye(dim)"
              1          1          1

##  g = 1  - circulant                                   'circulant'
              1          1           1

## g = -1  - circulant                                   'anticirculant'
              1          0           1

## g-Circulant                                             ' g!=+1 & g!=-1 '
              0          0           0


_________________________________________________________




## Copyright (C) 2000 - 2008  Rolf Fabian
##
## This file is part of Octave.
##
## It is published under the 'GNU General Public License' of the 'Free
## Software Foundation', either Version 3 of the license, or any later
## version. This work is free software. You can redistribute it and/or
## modify it under the terms of the specified license.
##
## This work 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.
##
## In case you have not received a copy of the specified license along
## with this file, consult GNU-website <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} address@hidden = circulant
(@var{v},@var{g},@var{mode})
## Return the @address@hidden square address@hidden@var{n}] matrix
@var{x}
## constructed from vector @var{v} of length @var{n}.
##
## In default @var{mode}='c', @var{v} is contained in first column of
@var{x}.
## Subsequently, it is circular shifted by @var{g} elements in each step
while
## moving through increasing column indices. If @var{g} is a negative
integer,
## the applied shifts are anticircular. By default, @address@hidden
(circular)
## is used, which might by also invoked by empty @var{g} input.
##
## If @var{mode}='r' is specified, vector @var{v} appears along the first
row
## of @var{x} and its circular (or anticircular) @var{g}-shifted
counterparts
## are aligned along rows of increasing indices.
##
## The results for special cases @address@hidden and @address@hidden are
## usually called 'circulant' (instead of '1-circulant') and 'anticirculant'
## (instead of '-1-circulant'). They form special Toeplitz (@var{g}=+1) and
## special Hankel (@var{g}=-1) matrices.
##
## Any circulant as well as any real or pure imaginary anticirculant
matrices are
## known to be 'normal', i.e. they commute with their own complex conjugated
## transposed matrix. This is of importance because the eigenvectors of
'normal'
## matrices always form an orthogonal basis spanning complete
address@hidden@var{n}]
## vector space. Then @var{x} can be always unitarily diagonalized. 
@var{g}-
## circulant matrices @var{x} constructed from other @var{g}-values
generally
## do not pocess the 'normal' property.
## @seealso{shift, isnormal}
## @end deftypefn

## Author: Rolf Fabian <r dot fabian at jacobs-univerity dot de>
## last revision 2008-02-19


function x = circulant (v,g,mode)

if (nargin < 1 || nargin > 3)
   usage ("x = circulant (v, g, mode)");
endif

if ( isempty (v) )
  y = [];
  return;
endif

[r,c] = size (v = v(:).');
if ( all ([r,c] > 1) )
  error ("circulant: 1st argument needs to be a vector.");
endif

if (nargin < 3)
  mode = 'c';
  ## circulate along columns by default
endif

if ( ! ischar (mode) || length (mode) != 1 
  || ! ( mode == 'c' || mode == 'r' ) )
  error ("circulant: mode = 'r'|'c' required.");
endif

if (nargin < 2 || isempty (g))
  g = 1;
endif

if (! isscalar(g))
  error ("circulant: shift parameter g needs to be scalar +/- integer.");
endif

if (g == 1)
  ## catch default g case to avoid creation of
  ## unnecessary large ind matrices.
  ind = zeros (c,c);
  ## preallocate memory
  ind (1:(c-1),:) = reshape ( kron (ones (c-1,1),(1:c).'), c-1,c );
  ind ( c,: ) = c:-1:1;

else
  j = -g;

  if (j < -1)
      j .+= c;
  endif

  ind = reshape ( kron (ones (1,c+j), 1:c), c+j,c );
  ind ( (c+1):(c+j),: ) = [];
endif

x = reshape ( v(ind(:)), c,c );

if (mode == 'r')
  x = x.';
endif

endfunction


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

-- 
View this message in context: 
http://www.nabble.com/LinAlg%3A-matrix-normality-check-proposed%2C-g-circulant-matrices-tp15589808p15589808.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.