matrix SVD



Robert Dodier wrote:
> On 10/19/07, Ruben Henner Zilibowitz <rzilibowitz at yahoo.com.au> wrote:
>
>   
>> I was wondering if there's a built in maxima function to compute the
>> SVD of a matrix (ie the Singular Value Decomposition). Is there one?
>>     
>
> There is dgesvd in the lapack package which computes SVD
> for numerical arguments. There is no symbolic SVD.
>
>   
You can define your own SVD operator.  I have the following three.

load(linearalgebra)$

put(infolevel, false, linearalgebra)$

Sing_vals(A):= reverse(sort(sqrt(eigens_by_jacobi(transpose(A).A))[1]))$

Sing_vals_a(A):= 
reverse(sort(float(sqrt(realpart(eigenvalues(transpose(A).A))))[1]))$ 

Sym_sing_values(A):= 
reverse(sort(sqrt(realpart(eigenvalues(transpose(A).A)))[1]))$ 

alias(svd_a, Sing_vals_a)$

alias(svd,Sing_vals)$

alias(sym_svd,Sym_sing_values)$

Sing_vals_usg: prt("Sing_vals(A) gives the list of singular values of
the square matrix A using the eigens_by_jacobi algorithm.  This is the 
list of square-roots of the
eigenvalues of transpose(A).A
An alias is svd(A). See also svd_a which uses the standard eigenvalues
routine.  sym_svd(A) does it symbolically. ")$

Of these, "svd" works fairly fast.  The others are slow. 

-sen