function F = factorize (A,strategy,burble)
%FACTORIZE an object-oriented method for solving linear systems
% and least-squares problems, and for representing operations with the
% inverse of a square matrix or the pseudo-inverse of a rectangular matrix.
%
% F = factorize(A) returns an object F that holds the factorization of A.
% x=F\b then solves a linear system or a least-squares problem. S=inverse(F)
% or S=inverse(A) returns a factorized representation of the inverse of A so
% that inverse(A)*b is mathematically equivalent to pinv(A)*b, but the former
% does not actually compute the inverse or pseudo-inverse of A.
%
% Example
%
% F = factorize(A) ; % LU, QR, or Cholesky factorization of A
% x = F\b ; % solve A*x=b; same as x=A\b
% S = inverse (F) ; % S represents the factorization of inv(A)
% x = S*b ; % same as x = A\b.
% E = A-B*inverse(D)*C % efficiently computes the Schur complement
% E = A-B*inv(D)*C % bad method for computing the Schur complement
% S = inverse(A) ; S(:,1) % compute just the first column of inv(A),
% % without computing inv(A)
%
% F = factorize (A, strategy, burble) ; % optional 2nd and 3rd inputs
%
% A string can be specified as a second input parameter to select the strategy
% used to factorize the matrix. The first two are meta-strategies:
%
% 'default' if rectangular
% use QR for A or A' (whichever is tall and thin).
% else
% if symmetric
% if positive real diagonal: try CHOL
% else (or if CHOL fails): try LDL
% end
% if not yet factorized: try LU (fails if rank-deficient)
% end
% if all else fails, or if QR or LU report that the matrix
% is singular (or nearly so): use COD
% This strategy mimics backslash, except that backslash never
% uses COD. Backslash also exploits other solvers, such as
% specialized tridiagonal and banded solvers.
%
% 'symmetric' as 'default' above, but assumes A is symmetric without
% checking, which is faster if you already know A is symmetric.
% Uses tril(A) and assumes triu(A) is its transpose. Results
% will be incorrect if A is not symmetric. If A is rectangular,
% the 'default' strategy is used instead.
%
% 'unsymmetric' as 'default', but assumes A is unsymmetric.
%
% The next "strategies" just select a single method, listed in decreasing order
% of generality and increasing order of speed and memory efficiency. All of
% them except the SVD can exploit sparsity.
%
% 'svd' use SVD. Never fails ... unless it runs out of time or memory.
% Coerces a sparse matrix A to full.
%
% 'cod' use COD. Almost as accurate as SVD, and much faster. Based
% on dense or sparse QR with rank estimation. Handles
% rank-deficient matrices, as long as it correctly estimates
% the rank. If the rank is ill-defined, use the SVD instead.
% Sparse COD requires the SPQR package to be installed
% (see http://www.cise.ufl.edu/research/sparse/SuiteSparse).
%
% 'qr' use QR. Reports a warning if A is singular.
%
% 'lu' use LU. Fails if A is rectangular; warning if A singular.
%
% 'ldl' use LDL. Fails if A is rank-deficient or not symmetric. Uses
% tril(A) and assumes triu(A) is the transpose of tril(A).
%
% 'chol' use CHOL. Fails if A is rank-deficient or not symmetric
% positive definite. If A is sparse, it uses tril(A) and
% assumes triu(A) is the transpose of tril(A). If A is dense,
% triu(A) is used instead.
%
% A third input, burble, can be provided to tell this function to print what
% methods it tries (if burble is nonzero).
%
% For a demo type "fdemo" in the Demo directory or see the Doc/ directory.
%
% See also inverse, slash, linsolve, spqr.
% Copyright 2011, Timothy A. Davis, University of Florida.
assert (ndims (A) == 2, 'Matrix must be 2D.') ;
if (nargin < 2 || isempty (strategy))
strategy = 'default' ;
end
if (nargin < 3)
burble = 0 ;
end
if (burble)
fprintf ('\nfactorize: strategy %s, A has size %d-by-%d, ', ...
strategy, size (A)) ;
if (issparse (A))
fprintf ('sparse with %d nonzeros.\n', nnz (A)) ;
else
fprintf ('full.\n') ;
end
end
switch strategy
case 'default'
[F me] = backslash_mimic (A, burble, 0) ;
case 'symmetric'
[F me] = backslash_mimic (A, burble, 1) ;
case 'unsymmetric'
[F me] = backslash_mimic (A, burble, 2) ;
case 'svd'
[F me] = factorize_svd (A, burble) ;
case 'cod'
[F me] = factorize_cod (A, burble) ;
case 'qr'
[F me] = factorize_qr (A, burble, 0) ;
case 'lu'
% do not report a failure if the matrix is singular
[F me] = factorize_lu (A, burble, 0) ;
case 'ldl'
[F me] = factorize_ldl (A, burble) ;
case 'chol'
[F me] = factorize_chol (A, burble) ;
otherwise
error ('FACTORIZE:invalidStrategy', 'unrecognized strategy.') ;
end
if (~isobject (F))
throw (me) ;
end
%-------------------------------------------------------------------------------
function [F me] = backslash_mimic (A, burble, strategy)
%BACKSLASH_MIMIC automatically select a method to factorize A.
F = [ ] ;
me = [ ] ;
[m n] = size (A) ;
% If the following condition is true, then the QR, QRT, or LU factorizations
% will report a failure if A is singular (or nearly so). This allows COD
% or COD_SPARSE to then be used instead. COD_SPARSE relies on the SPQR
% mexFunction in SuiteSparse, which might not be installed. In this case,
% QR, QRT, and LU do not report failures for sparse matrices that are singular
% (or nearly so), since there is no COD_SPARSE to fall back on.
fail_if_singular = ~issparse (A) || (exist ('spqr') == 3) ; %#ok
if (m ~= n)
% use QR for the rectangular case (ignore 'strategy' argument)
[F me] = factorize_qr (A, burble, fail_if_singular) ;
else
% square case: Cholesky, LDL, or LU factorization of A
switch strategy
case 0
is_symmetric = (nnz (A-A') == 0) ;
case 1
is_symmetric = true ;
case 2
is_symmetric = false ;
end
if (is_symmetric)
% A is symmetric (or assumed to be so)
d = diag (A) ;
if (all (d > 0) && nnz (imag (d)) == 0)
% try a Cholesky factorization
[F me] = factorize_chol (A, burble) ;
end
if (~isobject (F) && (~issparse (A) || isreal (A)))
% try an LDL factorization.
% complex sparse LDL does not yet exist in MATLAB
[F me] = factorize_ldl (A, burble) ;
end
end
if (~isobject (F))
% use LU if Cholesky and/or LDL failed, or were skipped.
[F me] = factorize_lu (A, burble, fail_if_singular) ;
end
end
if (~isobject (F))
% everything else failed, matrix is rank-deficient. Use COD
[F me] = factorize_cod (A, burble) ;
end
%-------------------------------------------------------------------------------
function [F me] = factorize_qr (A, burble, fail_if_singular)
% QR fails if the matrix is rank-deficient.
F = [ ] ;
me = [ ] ;
try
[m n] = size (A) ;
if (m >= n)
if (burble)
fprintf ('factorize: try QR of A ... ') ;
end
if (issparse (A))
F = factorization_qr_sparse (A, fail_if_singular) ;
else
F = factorization_qr_dense (A, fail_if_singular) ;
end
else
if (burble)
fprintf ('factorize: try QR of A'' ... ') ;
end
if (issparse (A))
F = factorization_qrt_sparse (A, fail_if_singular) ;
else
F = factorization_qrt_dense (A, fail_if_singular) ;
end
end
if (burble)
fprintf ('OK.\n') ;
end
catch me
if (burble)
fprintf ('failed.\nfactorize: %s\n', me.message) ;
end
end
%-------------------------------------------------------------------------------
function [F me] = factorize_chol (A, burble)
% LDL fails if the matrix is rectangular, rank-deficient, or not positive
% definite. Only the lower triangular part of A is used.
F = [ ] ;
me = [ ] ;
try
if (burble)
fprintf ('factorize: try CHOL ... ') ;
end
if (issparse (A))
F = factorization_chol_sparse (A) ;
else
F = factorization_chol_dense (A) ;
end
if (burble)
fprintf ('OK.\n') ;
end
catch me
if (burble)
fprintf ('failed.\nfactorize: %s\n', me.message) ;
end
end
%-------------------------------------------------------------------------------
function [F me] = factorize_ldl (A, burble)
% LDL fails if the matrix is rectangular or rank-deficient.
% Only the lower triangular part of A is used.
F = [ ] ;
me = [ ] ;
try
if (burble)
fprintf ('factorize: try LDL ... ') ;
end
if (issparse (A))
F = factorization_ldl_sparse (A) ;
else
F = factorization_ldl_dense (A) ;
end
if (burble)
fprintf ('OK.\n') ;
end
catch me
if (burble)
fprintf ('failed.\nfactorize: %s\n', me.message) ;
end
end
%-------------------------------------------------------------------------------
function [F me] = factorize_lu (A, burble, fail_if_singular)
% LU fails if the matrix is rectangular or rank-deficient.
F = [ ] ;
me = [ ] ;
try
if (burble)
fprintf ('factorize: try LU ... ') ;
end
if (issparse (A))
F = factorization_lu_sparse (A, fail_if_singular) ;
else
F = factorization_lu_dense (A, fail_if_singular) ;
end
if (burble)
fprintf ('OK.\n') ;
end
catch me
if (burble)
fprintf ('failed.\nfactorize: %s\n', me.message) ;
end
end
%-------------------------------------------------------------------------------
function [F me] = factorize_cod (A, burble)
% COD only fails when it runs out of memory.
F = [ ] ;
me = [ ] ;
try
if (burble)
fprintf ('factorize: try COD ... ') ;
end
if (issparse (A))
F = factorization_cod_sparse (A) ;
else
F = factorization_cod_dense (A) ;
end
if (burble)
fprintf ('OK.\n') ;
end
catch me
if (burble)
fprintf ('failed.\nfactorize: %s\n', me.message) ;
end
end
%-------------------------------------------------------------------------------
function [F me] = factorize_svd (A, burble)
% SVD only fails when it runs out of memory.
F = [ ] ;
me = [ ] ;
try
if (burble)
fprintf ('factorize: try SVD ... ') ;
end
F = factorization_svd (A) ;
if (burble)
fprintf ('OK.\n') ;
end
catch me
if (burble)
fprintf ('failed.\nfactorize: %s\n', me.message) ;
end
end
