functions in linalg.i -

             local sigma, u, vt; _sv_get_dcmp(a);  
 
       -or- local sigma, u, vt; _sv_get_dcmp(a, full)  
     Private routine used to extract Singular Value Decomposition of A into  
     external symbols: SIGMA, U and VT.  
       

             cholesky(a)  
 
       -or- cholesky(a, 0/1)  
         
     Given  a  symmetric  positive  definite  matrix  A,  returns  a  lower  
     triangular  matrix C  which is  the Cholesky  decomposition of  A such  
     that:  
       A = C(+,)*C(+,);  
     If  optional  second argument  is  true  (non-nil  and non-zero),  the  
     scratch values in  the upper triangular part of  C are left unchanged;  
     otherwise (the  default behavior), the  upper triangular part of  C is  
     filled with zeros.  

             diag(a)  
 
     Returns the diagonal of matrix A (if A is a 2D array) or a square  
     diagonal matrix with diagonal elements equal to those of vector A (if  
     A is a 1D array).  
       

             euclidean_norm(x)  
 
    Returns the Euclidian norm of X: sqrt(sum(X*X)), taking care of overflows.  

             gram_schmidt_orthonormalization, b;  
 
       -or- gram_schmidt_orthonormalization(b)  
     Performs Gram-Schmidt orthonormalization of basis functions B.  If  
     B is an array of pointers, then the input basis vectors are *B(i)  
     for i=1,..,numberof(B); otherwise, the input basis vectors are B(..,i)  
     for i=1,..,dimsof(B)(0).  When called as a subroutine, the operation  
     is done "in-place".  
       

             pm, a;  
 
     Prints matrix A.  Keyword FILE can be set to the name/stream of output  
     file; the default is to use standard output.  If keyword NDIGITS is  
     set, floating-point values get printed with that number of significant  
     digits; alternatively, keyword FORMAT may be set with the format for  
     each element of A -- for complex numbers, the real and imaginary parts  
     use the same format.  
       

             sv_dcmp(a)  
 
       -or- sv_dcmp(a, full)  
     Computes the singular value decomposition of matrix A.  The result  
     is an array of pointers:  
         [&SIGMA, &U, &VT]  
     where SIGMA is the vector of  singular values of A, and the columns of  
     U  and  the  rows of  VT  are  the  left  and right  singular  vectors  
     (respectively).  Using matrix notation:  
         A = U.diag(SIGMA).V' = U.diag(SIGMA).VT  
     and using Yorick notation:  
         A = (U*SIGMA(-,))(,+)*VT(+,)  
           = U(,+)*(SIGMA*VT)(+,)  
     FULL has the same meaning as keyword FULL in SVdec (to see).  
       

             sv_intro - introduction to SVD Yorick package  
 
    
  Notes about matrix multiplication in Yorick:  
    
    A.B  = A(,+)*B(+,)  
    A'.B = A(+,)*B(+,)  // transpose of A times B  
    A.B' = A(,+)*B(,+)  // A times transpose of B  
    
    diag(s).A  = diag(s)(,+)*A(+,)  
               = s*A = A*s  
    
    diag(s).A' = diag(s)(,+)*A(,+)  
               = transpose(s(-,)*A) = transpose(A*s(-,))  
    
    A.diag(s)  = A(,+)*diag(s)(+,)  
               = A*s(-,) = s(-,)*A  
    
    A'.diag(s) = A(+,)*diag(s)(+,) = A(+,)*diag(s)(,+)  
               = transpose(A*s) = transpose(s*A)  
    
  Singular Value Decomposition:  
    
    A = U.diag(SIGMA).V' = U.diag(SIGMA).VT  
      = (U*SIGMA(-,))(,+)*VT(+,)  
      = U(,+)*(SIGMA*VT)(+,)  
    
  where:  
    
    SIGMA = SVdec(A, U, VT)  
    
  Columns  of  U  and  V form  an  orthonormal basis  (i.e.  U  and V  are  
  column-orthonormal):  
    
    U'.U = V'.V = I  
    
  in Yorick notation:  
    
    U(+,)*U(+,) = V(+,)*V(+,) = VT(,+)*VT(,+) = unit(N)  
    
  Note (to be verified): if U and/or V are square, they are also  
  row-orthonormal:  
    
    U'.U = U.U' = I   (if U is square)  
    V'.V = V.V' = I   (if V is square)  
    
  Generalized-inverse of A:  
    
    AP = V.diag(1/SIGMA).U' = VT'.diag(1/SIGMA).U'  
      = VT(+,)*((1.0/SIGMA)(-,)*U)(,+)  
      = ((1.0/SIGMA)*VT)(+,)*U(,+)  
    
  Solving a linear problem: A.x ~ b with SVD (taking care of index  
  ordering for faster matrix multiplication):  
    
    X = V.diag(W).U'.B  
      = (W*VT)(+,)*U(,+))(,+)*B(+,..)  
      = (W*VT)(+,)*(U(+,)*B(+,..))(+,..)  // sum over 1st indices: faster  
    
  where W is an approximation of 1/SIGMA.  
    

             sv_solve_trunc(a, b, eps)  
 
     Solve linear problem A.x = b by truncated singular value method.  A is  
     either a matrix or the  singular value value decomposition as returned  
     by sv_dcmp  (to see).  B  is the right  hand side vector (or  array to  
     solve for  several right  hand sides  at a time).   EPS (in  the range  
     [0,1])  is the  minimum  relative  singular value  to  keep, i.e.  all  
     singular values  less than EPS*max(SIGMA) get discarded  (SIGMA is the  
     vector of singular values of A).  The result is:  
       
         (W*VT)(+,)*(U(+,)*B(+,..))(+,..)  
         
     where  SIGMA,  U and  VT  are the  components  of  the singular  value  
     decomposition of A and W is:  
         W(i) = 1/SIGMA(i)   if SIGMA(i) > EPS*max(SIGMA);  
                0            otherwise.  
       

             sv_solve_wiener(a, b, eps)  
 
     Solve  linear problem  A.x =  b by  Wiener filtering  of  the singular  
     values.   A   is  either  a   matrix  or  the  singular   value  value  
     decomposition as  returned by sv_dcmp (to  see).  B is  the right hand  
     side  vector (or  array to  solve for  several right  hand sides  at a  
     time).  EPS (in the range [0,1]) is the relative singular value filter  
     level.  The result is:  
       
         (W*VT)(+,)*(U(+,)*B(+,..))(+,..)  
         
     where  SIGMA,  U and  VT  are the  components  of  the singular  value  
     decomposition of A and W is:  
         W = SIGMA/(SIGMA^2 + (EPS*max(SIGMA))^2)  
       

             trace(a)  
 
     Returns the trace of matrix A.