sparseinv

SPARSEINV computes the sparse inverse subset of a real sparse square matrix A.
 This function is typically much faster than computing all of inv(A).

 [Z, Zpattern, L, D, U, P, Q, stats] = sparseinv (A)

 Z is a subset of the inverse a sparse matrix A of full rank.  On output, if
 Zpattern(i,j)=1, it means that Z(i,j) has been computed.  That is, the norm
 of (Zpattern .* (Z - inv (A))) will be small.

 Method: The permuted matrix C = P*A*Q is first factorized into C =
 (L+I)*D*(U+I) where D is diagonal, L+I is lower triangular with unit
 diagonal, and U+I is upper triangular with unit diagonal (I = speye (n)).  If
 A is symmetric and positive definite, then a Cholesky factorization is used
 (in which case P=Q' and L=U', and Z will include all diagonal entries of
 inv(A)).  Next, the entries in the inverse of C that correspond to nonzero
 values in Zpattern are found via Takahashi's method.  Zpattern is the
 symbolic Cholesky factorization of C+C', so it includes all entries in L+U
 and its transpose.

 stats is an optional struct containing statistics on the factorization.

 Example:
   load west0479
   A = west0479 ;
   [Z, Zpattern] = sparseinv (A) ;
   S = inv (A) ;
   err = norm (Zpattern .* (Z - S), 1) / norm (S, 1)

 See also inv, lu, chol.

0001 function [Z, Zpattern, L, D, U, P, Q, stats] = sparseinv (A)
0002 %SPARSEINV computes the sparse inverse subset of a real sparse square matrix A.
0003 % This function is typically much faster than computing all of inv(A).
0004 %
0005 % [Z, Zpattern, L, D, U, P, Q, stats] = sparseinv (A)
0006 %
0007 % Z is a subset of the inverse a sparse matrix A of full rank.  On output, if
0008 % Zpattern(i,j)=1, it means that Z(i,j) has been computed.  That is, the norm
0009 % of (Zpattern .* (Z - inv (A))) will be small.
0010 %
0011 % Method: The permuted matrix C = P*A*Q is first factorized into C =
0012 % (L+I)*D*(U+I) where D is diagonal, L+I is lower triangular with unit
0013 % diagonal, and U+I is upper triangular with unit diagonal (I = speye (n)).  If
0014 % A is symmetric and positive definite, then a Cholesky factorization is used
0015 % (in which case P=Q' and L=U', and Z will include all diagonal entries of
0016 % inv(A)).  Next, the entries in the inverse of C that correspond to nonzero
0017 % values in Zpattern are found via Takahashi's method.  Zpattern is the
0018 % symbolic Cholesky factorization of C+C', so it includes all entries in L+U
0019 % and its transpose.
0020 %
0021 % stats is an optional struct containing statistics on the factorization.
0022 %
0023 % Example:
0024 %   load west0479
0025 %   A = west0479 ;
0026 %   [Z, Zpattern] = sparseinv (A) ;
0027 %   S = inv (A) ;
0028 %   err = norm (Zpattern .* (Z - S), 1) / norm (S, 1)
0029 %
0030 % See also inv, lu, chol.
0031 
0032 % Copyright 2011, Timothy A. Davis, http://www.suitesparse.com
0033 
0034 get_stats = (nargout > 7) ;
0035 if (get_stats)
0036     t1 = tic ;
0037 end
0038 
0039 % check inputs
0040 if (~issparse (A))
0041     error ('A must be sparse') ;
0042 end
0043 [m n] = size (A) ;
0044 if (m ~= n)
0045     error ('A must be square') ;
0046 end
0047 if (~isreal (A))
0048     error ('complex matrices not supported') ;
0049 end
0050 
0051 % construct the factorization: C = P*A*Q = (L+I)*D*(U+I)
0052 p = 1 ;
0053 if (all (diag (A)) > 0 && nnz (A-A') == 0)
0054     [L,p,Q] = chol (A, 'lower') ;
0055 end
0056 if (p == 0)
0057     % Cholesky worked.
0058     P = Q' ;
0059     d = diag (L) ;
0060     L = tril (L / diag (d), -1) ;
0061     U = L' ;
0062     d = d.^2 ;
0063     D = diag (d) ;
0064 else
0065     % Cholesky failed, or wasn't attempted.  Use LU instead.
0066     [L,U,P,Q] = lu (A) ;
0067     d = diag (U) ;
0068     if (any (d == 0))
0069         error ('A must be full-rank') ;
0070     end
0071     D = diag (d) ;
0072     U = triu (D \ U, 1) ;
0073     L = tril (L, -1) ;
0074 end
0075 d = full (d) ;
0076 
0077 % find the symbolic Cholesky of C+C'
0078 S = spones (P*A*Q) ;
0079 [c,h,pa,po,R] = symbfact (S+S') ;
0080 clear h pa po
0081 Zpattern = spones (R+R') ;
0082 clear R S
0083 
0084 if (get_stats)
0085     t1 = toc (t1) ;
0086     t2 = tic ;
0087 end
0088 
0089 % compute the sparse inverse subset
0090 [Z takflops] = sparseinv_mex (L, d, U', Zpattern) ;
0091 if (p == 0)
0092     % Force Z to be symmetric.  This is because sparseinv_mex does not
0093     % exploit the symmetry in the factorization, but computes both upper and
0094     % lower triangular parts of Z separately.  The work for the Takahashi
0095     % equations could be cut in half as a result.
0096     Z = (Z + Z') / 2 ;
0097 end
0098 
0099 % permute the result
0100 Z = Q*Z*P ;
0101 Zpattern = Q*Zpattern*P ;
0102 
0103 % return stats, if requested
0104 if (nargout > 7)
0105     t2 = toc (t2) ;
0106     if (p == 0)
0107         stats.kind = 'Cholesky' ;
0108         fl = 2*n + sum (c.^2) ;
0109         stats.nnz_factors = sum (c) ;
0110     else
0111         stats.kind = 'LU' ;
0112         Lnz = full (sum (spones (L))) ;            % off diagonal nz in cols of L
0113         Unz = full (sum (spones (U')))' ;    % off diagonal nz in rows of U
0114         fl = n + 2*Lnz*Unz + sum (Lnz) ;
0115         stats.nnz_factors = nnz (L) + nnz (U) + n ;
0116     end
0117     stats.flops_factorization = fl ;
0118     stats.flops_Takahashi = takflops ;
0119     stats.time_factorization = t1 ;
0120     stats.time_Takahashi = t2 ;
0121 end
0122