Re: hess(rosser())

[Ed Meyer]

On Wed, Jul 4, 2012 at 10:30 AM, Matyas Sustik <address@hidden> wrote:

On Wed, Jul 4, 2012 at 12:12 PM, Juan Pablo Carbajal <address@hidden> wrote:

I did some minimal work on manipulation of symmetric matrices (not
optimizing storage though).

- There are the functions vech and unvech (OF package general-1.3.1)
are for converting from (anti)symmetric to vector form
http://octave.svn.sf.net/viewvc/octave/trunk/octave-forge/main/general/inst/unvech.m?revision=10727&view=markup

- The function unvech has a subfunction sub2ind_tril (this one and
ind2sub_tril are also private in the OF package mechanics
http://octave.svn.sf.net/viewvc/octave/trunk/octave-forge/main/mechanics/inst/core/private/).
JordiGh gave me some ideas on how to implement this in the source code
of sub2ind and ind2sub functions, but I haven't done it yet. A short
explanation can be found in my TODO list in the wiki
http://octave.org/wiki/index.php?title=User:KaKiLa#extend_sub2ind_and_ind2sub

My two cents.

In my effort I would restrict myself to making sure that the most suited LAPACK function is being called when a matrix is known symmetric or Hermitian.  That is:

hess() would call dsytrd, then

eig() would call dstedc (maybe dsteqr or dstemr depending on parm?) on the resulting tridiagonal matrix.

(For Hermitian hess() would call chetrd etc.)

What is not clear to me is how octave keeps track whether a matrix is symmetric or Hermitian.  For example rosser() returns a symmetric matrix, but I am not sure that it is "marked" as symmetric.  Checking every time that A= A' would be inefficient.

This improvement stems from the observation that hess(rosser()) returns a matrix that is (quite) not symmetric while rosser() is.

-Matyas


 

 

--
M. Sc. Juan Pablo Carbajal
-----
PhD Student
University of Zürich
http://ailab.ifi.uzh.ch/carbajal/

>> I looked at the source code.  It appears that (dense) symmetric matrices are not
>> handled specially in the eigenvalue routines.  However, the concept of symmetry in >> octave does exist.

Looking at the source it appears to me (verified with a debugger) that symmetry
is accounted for in eig() and the hessenberg routines are not called except from
hess(). Since a Hessenberg reduction is intended for nonsymmetric matrices
(though it will produce a tridiagonal matrix in the symmetric case with exact arithmetic)
it seems like an unnecessary complication to make hess() do something it was not
intended for and currently it shows what the Hessenberg reduction does to a symmetric matrix with non-exact arithmetic.


-- 
Ed Meyer
206-282-3160 (H)
206-453-8778 (M)