Code covered by the BSD License

# Don't let that INV go past your eyes; to solve that system, FACTORIZE!

### Tim Davis (view profile)

14 May 2009 (Updated )

A simple-to-use object-oriented method for solving linear systems and least-squares problems.

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

factorization
```classdef factorization
%FACTORIZATION a generic matrix factorization object
% Normally, this object is created via the F=factorize(A) function.  Users
% do not need to use this method directly.
%
% This is an abstract class that is specialized into 13 different kinds of
% matrix factorizations:
%
%   factorization_chol_dense    dense Cholesky      A = R'*R
%   factorization_lu_dense      dense LU            A(p,:) = L*U
%   factorization_qr_dense      dense QR of A       A = Q*R
%   factorization_qrt_dense     dense QR of A'      A' = Q*R
%   factorization_ldl_dense     dense LDL           A(p,p) = L*D*L'
%   factorization_cod_dense     dense COD           A = U*R*V'
%
%   factorization_chol_sparse   sparse Cholesky     P'*A*P = L*L'
%   factorization_lu_sparse     sparse LU           P*(R\A)*Q = L*U
%   factorization_qr_sparse     sparse QR of A      (A*P)'*(A*P) = R'*R
%   factorization_qrt_sparse    sparse QR of A'     (P*A)*(P*A)' = R'*R
%   factorization_ldl_sparse    sparse LDL          P'*A*P = L*D*L'
%   factorization_cod_sparse    sparse COD          A = U*R*V'
%
%   factorization_svd           SVD                 A = U*S*V'
%
% The abstract class provides the following functions.  In the descriptions,
% F is a factorization.  The arguments b, y, and z may be factorizations or
% matrices.  The output x is normally matrix unless it can be represented as a
% scaled factorization.  For example, G=F\2 and G=inverse(F)*2 both return a
% factorization G.  Below, s is always a scalar, and C is always a matrix.
%
%   These methods return a matrix x, unless one argument is a scalar (in which
%   case they return a scaled factorization object):
%   x = mldivide (F, b)     x = A \ b
%   x = mrdivide (b, F)     x = b / A
%   x = mtimes (y, z)       y * z
%
%   These methods always return a factorization:
%   F = uplus (F)           +F
%   F = uminus (F)          -F
%   F = inverse (F)         representation of inv(A), without computing it
%   F = ctranspose (F)      F'
%
%   These built-in methods return a scalar:
%   s = isreal (F)
%   s = isempty (F)
%   s = isscalar (F)
%   s = issingle (F)
%   s = isnumeric (F)
%   s = isfloat (F)
%   s = isvector (F)
%   s = issparse (F)
%   s = isfield (F,f)
%   s = isa (F, s)
%   s = condest (F)
%
%   This method returns the estimated rank from the factorization.
%   s = rankest (F)
%
%   These methods support access to the contents of a factorization object
%   e = end (F, k, n)
%   [m,n] = size (F, k)
%   S = double (F)
%   C = subsref (F, ij)
%   S = struct (F)
%   disp (F)
%
% The factorization_chol_dense object also provides cholupdate, which acts
% just like the builtin cholupdate.
%
% The factorization_svd object provides:
%
%   c = cond (F,p)      the p-norm condition number.  p=2 is the default.
%                       cond(F,2) takes no time to compute, since it was
%                       computed when the SVD factorization was found.
%   a = norm (F,p)      the p-norm.  see the cond(F,p) discussion above.
%   r = rank (F)        returns the rank of A, precomputed by the SVD.
%   Z = null (F)        orthonormal basis for the null space of A
%   Q = orth (F)        orthonormal basis for the range of A
%   C = pinv (F)        the pseudo-inverse, V'*(S\V).
%   [U,S,V] = svd (F)   SVD of A or pinv(A), regular, economy, or rank-sized
%

% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com

properties (SetAccess = protected)
% The abstract class holds a QR, LU, Cholesky factorization:
A = [ ] ;           % a copy of the input matrix
Factors = [ ] ;
is_inverse = false ;% F represents the factorization of A or inv(A)
is_ctrans = false ; % F represents the factorization of A or A'
kind = '' ;         % a string stating what kind of factorization F is
alpha = 1 ;         % F represents the factorization of A or alpha*A.
A_rank = [ ] ;      % rank of A, from SVD, or estimate from COD
A_cond = [ ] ;      % 2-norm condition number of A, from SVD
A_condest = [ ] ;   % quick and dirty estimate of the condition number
% If F is inverted, alpha doesn't change.  For example:
%   F = alpha*factorize(A) ; % F = alpha*A, in factorized form.
%   G = inverse(F) ;         % G = inv(alpha*A)
%   H = beta*G               % H = beta*inv(alpha*A)
%                                = inv((alpha/beta)*A)
% So to update alpha via scaling, beta*F, the new scale factor beta
% is multiplied with F.alpha if F.is_inverse is false.  Otherwise,
% F.alpha is divided by beta to get the new scale factor.
end

methods (Abstract)
x = mldivide_subclass (F, b) ;
x = mrdivide_subclass (b, F) ;
e = error_check (F) ;
end

methods

%-----------------------------------------------------------------------
% mldivide and mrdivide: return a scaled factorization or a matrix
%-----------------------------------------------------------------------

% Let b be a double scalar, F a non-scalar factorization, and g a scalar
% factorization.  Then these operations return scaled factorization
% objects (unless flatten is true, in which case a matrix is returned):
%
%   F\b = inverse (F) * b           |   F/b = F / b
%   F\g = inverse (F) * double (g)  |   F/g = F / double (g)
%   b\F = F / b                     |   b/F = b * inverse (F)
%   g\F = F / double (g)            |   g/F = double (g) * inverse (F)
%
% Otherwise mldivide & mrdivide always return a matrix as their result.

function x = mldivide (y, z, flatten)
%MLDIVIDE x = y\z where either y or z or both are factorizations.
flatten = (nargin > 2 && flatten) ;
if (isobject (y) && isscalar (z) && ~flatten)
% x = y\z where y is an object and z is scalar (perhaps object).
% result is a scaled factorization object x.
x = scale_factor (inverse (y), ~(y.is_inverse), double (z)) ;
elseif (isscalar (y) && isobject (z) && ~flatten)
% x = y\z where y is scalar (perhaps object) and z is an object.
% result is a scaled factorization object x.
x = scale_factor (z, ~(z.is_inverse), double (y)) ;
else
% result x will be a matrix.  b is coerced to be a matrix.
[F, b, first_arg_is_F] = getargs (y, z) ;
if (~first_arg_is_F)
% x = b\F where F is a factorization.
error_if_inverse (F, 1) ;
x = b \ F.A ;            % use builtin backslash
elseif (F.is_ctrans)
% x = F\b where F represents (alpha*A)' or inv(alpha*A)'
if (F.is_inverse)
% x = inv(alpha*A)'\b = (A'*b)*alpha'
x = (F.A'*b) * F.alpha' ;
else
% x = (alpha*A)'\b = (b'/A)' / alpha'
x = mrdivide_subclass (b', F)' / F.alpha' ;
end
else
% x = F\b where F represents (alpha*A) or inv(alpha*A)
if (F.is_inverse)
% x = inv(alpha*A)\b = (A*b)*alpha
x = (F.A*b) * F.alpha ;
else
% x = (alpha*A)\b = (A\b) / alpha
x = mldivide_subclass (F, b) / F.alpha ;
end
end
end
end

function x = mrdivide (y, z, flatten)
%MRDIVIDE x = y/z where either y or z or both are factorizations.
flatten = (nargin > 2 && flatten) ;
if (isobject (y) && isscalar (z) && ~flatten)
% x = y/z where y is an object and z is scalar (perhaps object).
% result is a scaled factorization object x.
x = scale_factor (y, ~(y.is_inverse), double (z)) ;
elseif (isscalar (y) && isobject (z) && ~flatten)
% x = y/z where y is scalar (perhaps object) and z is an object.
% result is a scaled factorization object x.
x = scale_factor (inverse (z), ~(z.is_inverse), double (y)) ;
else
% result x will be a matrix.  b is coerced to be a matrix.
[F, b, first_arg_is_F] = getargs (y, z) ;
if (first_arg_is_F)
% x = F/b where F is a factorization object
error_if_inverse (F, 2)
x = F.A / b ;            % use builtin slash
elseif (F.is_ctrans)
% x = b/F where F represents (alpha*A)' or inv(alpha*A)'
if (F.is_inverse)
% x = b/inv(alpha*A)' = (b*A')*alpha'
x = (b*F.A') * F.alpha' ;
else
% x = b/(alpha*A)' = (A\b')' / alpha'
x = mldivide_subclass (F, b')' / F.alpha' ;
end
else
% x = b/F where F represents (alpha*A) or inv(alpha*A)
if (F.is_inverse)
% x = b/inv(alpha*A) = (b*A)*alpha
x = (b*F.A) * F.alpha ;
else
% x = b/(alpha*A) = (b/A) / alpha
x = mrdivide_subclass (b, F) / F.alpha ;
end
end
end
end

%-----------------------------------------------------------------------
% mtimes: a simple and clean wrapper for mldivide and mrdivide
%-----------------------------------------------------------------------

function x = mtimes (y, z)
%MTIMES x=y*z where y or z is a factorization object (or both).
% Since inverse(F) is so cheap, and does the right thing inside
% mldivide and mrdivide, this is just a simple wrapper.
if (isobject (y))
% A*b               becomes inverse(A)\b
% inverse(A)*b      becomes A\b
% A'*b              becomes inverse(A)'\b
% inverse(A)'*b     becomes A'\b
x = mldivide (inverse (y), z) ;
else
% b*A               becomes b/inverse(A)
% b*inverse(A)      becomes b/A
% b*A'              becomes b/inverse(A)'
% b*inverse(A)'     becomes b/A'
% y is a scalar or matrix, z must be an object
x = mrdivide (y, inverse (z)) ;
end
end

%-----------------------------------------------------------------------
% uplus, uminus, ctranspose, inverse: always return a factorization
%-----------------------------------------------------------------------

function F = uplus (F)
%UPLUS +F
end

function F = uminus (F)
%UMINUS -F
F.alpha = -(F.alpha) ;
end

function F = inverse (F)
%INVERSE "inverts" F by flagging it as factorization of inv(A)
F.is_inverse = ~(F.is_inverse) ;
end

function F = ctranspose (F)
%CTRANSPOSE "transposes" F by flagging it as factorization of A'
F.is_ctrans = ~(F.is_ctrans) ;
end

%-----------------------------------------------------------------------
% is* methods that return a scalar
%-----------------------------------------------------------------------

function s = isreal (F)
%ISREAL for F=factorize(A): same as isreal(A)
s = isreal (F.A) ;
end

function s = isempty (F)
%ISEMPTY for F=factorize(A): same as isempty(A)
s = any (size (F.A) == 0) ;
end

function s = isscalar (F)
%ISSCALAR for F=factorize(A): same as isscalar(A)
s = isscalar (F.A)  ;
end

function s = issingle (F)                                           %#ok
%ISSINGLE for F=factorize(A) is always false
s = false ;
end

function s = isnumeric (F)                                          %#ok
%ISNUMERIC for F=factorize(A) is always true
s = true ;
end

function s = isfloat (F)                                            %#ok
%ISFLOAT for F=factorize(A) is always true
s = true ;
end

function s = isvector (F)
%ISVECTOR for F=factorize(A): same as isvector(A)
s = isvector (F.A) ;
end

function s = issparse (F)
%ISSPARSE for F=factorize(A): same as issparse(A)
s = issparse (F.A) ;
end

function s = isfield (F, f)                                         %#ok
%ISFIELD isfield(F,f) is true if F.f exists, false otherwise
s = (ischar (f) && (strcmp (f, 'A') ...
|| strcmp (f, 'Factors') || strcmp (f, 'kind') ...
|| strcmp (f, 'is_inverse') || strcmp (f, 'is_ctrans') ...
|| strcmp (f, 'alpha') || strcmp (f, 'A_rank') ...
|| strcmp (f, 'A_cond') || strcmp (f, 'A_condest'))) ;
end

function s = isa (F, s)
%ISA for F=factorize(A): 'double', 'numeric', 'float' are true.
% For other types, the builtin isa does the right thing.
s = strcmp (s, 'double') || strcmp (s, 'numeric') ||  ...
strcmp (s, 'float') || builtin ('isa', F, s) ;
end

%-----------------------------------------------------------------------
% condest, rankest
%-----------------------------------------------------------------------

function C = abs (F)
%ABS abs(F) returns abs(A) or abs(inverse(A)), as appropriate.  The
% ONLY reason abs is included here is to support the builtin
% normest1 for small matrices (n <= 4).  Computing abs(inverse(A))
% explicitly computes the inverse of A, so use with caution.
C = abs (double (F)) ;
end

function s = condest (F)
%CONDEST 1-norm condition number for square matrices.
% Does not require another factorization of A, so it's very fast.
% Does NOT explicitly compute the inverse of A.  Instead, if F
% represents an inverse, F*x inside normest1 does the right thing,
% and does A\b using the factorization F.
A = F.A ;                                                       %#ok
[m, n] = size (A) ;                                             %#ok
if (m ~= n)
error ('MATLAB:condest:NonSquareMatrix', ...
'Matrix must be square.') ;
end
if (n == 0)
s = 0 ;
elseif (F.is_inverse)
% F already represents the factorization of the inverse of A
s = F.alpha * norm (A,1) * normest1 (F) ;                   %#ok
else
% Note that the inverse is NOT explicitly computed.
s = F.alpha * norm (A,1) * normest1 (inverse (F)) ;         %#ok
end
end

function r = rankest (F)
%RANKEST returns the estimated rank of A.
% It is a very rough estimate if Cholesky, LU, QR, or LDL succeeded
% (in which A is assumed to have full rank).  COD finds a more
% accurate estimate, and SVD finds the exact rank.
r = F.A_rank ;
end

%-----------------------------------------------------------------------
% end, size
%-----------------------------------------------------------------------

function e = end (F, k, n)
%END returns index of last item for use in subsref
if (n == 1)
e = numel (F.A) ;   % # of elements, for linear indexing
else
e = size (F, k) ;   % # of rows or columns in A or pinv(A)
end
end

function [m, n] = size (F, k)
%SIZE returns the size of the matrix F.A in the factorization F
if (F.is_inverse ~= F.is_ctrans)
% swap the dimensions to match pinv(A)
if (nargout > 1)
[n, m] = size (F.A) ;
else
m = size (F.A) ;
m = m ([2 1]) ;
end
else
if (nargout > 1)
[m, n] = size (F.A) ;
else
m = size (F.A) ;
end
end
if (nargin > 1)
m = m (k) ;
end
end

%-----------------------------------------------------------------------
% double: a wrapper for subsref
%-----------------------------------------------------------------------

function S = double (F)
%DOUBLE returns the factorization as a matrix, A or inv(A)
ij.type = '()' ;
ij.subs = cell (1,0) ;
S = subsref (F, ij) ;   % let factorize.subsref do all the work
end

%-----------------------------------------------------------------------
% subsref: returns a matrix
%-----------------------------------------------------------------------

function C = subsref (F, ij)
%SUBSREF A(i,j) or (i,j)th entry of inv(A) if F is inverted.
% Otherwise, explicit entries in the inverse are computed.
% This method also extracts the contents of F with F.whatever.
switch (ij (1).type)
case '.'
% F.A usage: extract one of the matrices from F
switch ij (1).subs
case 'A'
C = F.A ;
case 'Factors'
C = F.Factors ;
case 'is_inverse'
C = F.is_inverse ;
case 'is_ctrans'
C = F.is_ctrans ;
case 'kind'
C = F.kind ;
case 'alpha'
C = F.alpha ;
case 'A_cond'
C = F.A_cond ;
case 'A_condest'
C = F.A_condest ;
case 'A_rank'
C = F.A_rank ;
otherwise
error ('MATLAB:nonExistentField', ...
'Reference to non-existent field ''%s''.', ...
ij (1).subs) ;
end
% F.X(2,3) usage, return X(2,3), for component F.X
if (length (ij) > 1 && ~isempty (ij (2).subs))
C = subsref (C, ij (2)) ;
end
case '()'
C = subsref_paren (F, ij) ;
case '{}'
error ('MATLAB:cellRefFromNonCell', ...
'Cell contents reference from a non-cell array object.') ;
end
end

%-----------------------------------------------------------------------
% struct: extracts all contents of a factorization object
%-----------------------------------------------------------------------

function S = struct (F)
%STRUCT convert factorization F into a struct.
% S cannot be used for subsequent object methods here.
S.A = F.A ;
S.Factors = F.Factors ;
S.is_inverse = F.is_inverse ;
S.is_ctrans = F.is_ctrans ;
S.alpha = F.alpha ;
S.A_rank = F.A_rank ;
S.A_cond = F.A_cond ;
S.kind = F.kind ;
end

%-----------------------------------------------------------------------
% disp: displays the contents of F
%-----------------------------------------------------------------------

function disp (F)
%DISP displays a factorization object
fprintf ('  class: %s\n', class (F)) ;
fprintf ('  %s\n', F.kind) ;
fprintf ('  A: [%dx%d double]\n', size (F.A)) ;
fprintf ('  Factors:\n') ; disp (F.Factors) ;
fprintf ('  is_inverse: %d\n', F.is_inverse) ;
fprintf ('  is_ctrans: %d\n', F.is_ctrans) ;
fprintf ('  alpha: %g', F.alpha) ;
if (~isreal (F.alpha))
fprintf (' + (%g)i', imag (F.alpha)) ;
end
fprintf ('\n') ;
if (~isempty (F.A_rank))
fprintf ('  A_rank: %d\n', F.A_rank) ;
end
if (~isempty (F.A_condest))
fprintf ('  A_condest: %d\n', F.A_condest) ;
end
if (~isempty (F.A_cond))
fprintf ('  A_cond: %d\n', F.A_cond) ;
end
end
end

%---------------------------------------------------------------------------
% methods that are not user-callable
%---------------------------------------------------------------------------

methods (Access = protected)

function [F, b, first_arg_is_F] = getargs (y, z)
first_arg_is_F = isobject (y) ;
if (first_arg_is_F)
F = y ;             % first argument is a factorization object
b = double (z) ;    % 2nd one coerced to be a matrix
else
b = y ;             % first argument is not an object
F = z ;             % second one must be an object
end
end

function F = scale_factor (F, use_beta_inverse, beta)
%SCALE_FACTOR scales a factorization
if (use_beta_inverse)
% F = inv(alpha*A), so F*beta = inv((alpha/beta)*A)
if (F.is_ctrans)
F.alpha = F.alpha / beta' ;
else
F.alpha = F.alpha / beta ;
end
else
% F = alpha*A, so F*beta = (alpha*beta)*A
if (F.is_ctrans)
F.alpha = F.alpha * beta' ;
else
F.alpha = F.alpha * beta ;
end
end
end
end
end

%-------------------------------------------------------------------------------
% subsref_paren: support function for subsref
%-------------------------------------------------------------------------------

function C = subsref_paren (F, ij)
%SUBSREF_PAREN C = subsref_paren(F,ij) implements C=F(i,j) and C=F(i)

% F(2,3) usage, return A(2,3) or the (2,3) of inv(A).
assert (length (ij) == 1, 'Improper index matrix reference.') ;
A = F.A ;
is_ctrans = F.is_ctrans ;
if (is_ctrans && length (ij.subs) > 1)   % swap i and j
ij.subs = ij.subs ([2 1]) ;
end

if (F.is_inverse)

% requesting explicit entries of the inverse

if (length (ij.subs) == 1)
% for linear indexing of the inverse (C=F(i)), first
% convert to double and then use builtin subsref
C = subsref (double (F), ij) ;
else
% standard indexing, C = F(i,j)
if (is_ctrans)
[n, m] = size (A) ;
else
[m, n] = size (A) ;
end
if (length (ij.subs) == 2)
ilen = length (ij.subs {1}) ;
if (strcmp (ij.subs {1}, ':'))
ilen = n ;
end
jlen = length (ij.subs {2}) ;
if (strcmp (ij.subs {2}, ':'))
jlen = m ;
end
j = ij ;
j.subs {1} = ':' ;
i = ij ;
i.subs {2} = ':' ;
if (jlen <= ilen)
% compute X=S(:,j) of S=inv(A) and return X(i,:)
C = subsref (mldivide (...
inverse (F), ...
subsref (identity (A, m), j), 1), i) ;
else
% compute X=S(i,:) of S=inv(A) and return X(:,j)
C = subsref (mrdivide (...
subsref (identity (A, n), i), ...
inverse (F), 1), j) ;
end
else
% the entire inverse has been explicitly computed
C = mldivide (inverse (F), identity (A, m), 1) ;
end
end

else

% F is not inverted, so just return A(i,j)
if (isempty (ij (1).subs))
C = A ;
else
C = subsref (A, ij) ;
end
C = C * F.alpha ;
if (is_ctrans)
C = C' ;
end
end
end

%-------------------------------------------------------------------------------
% identity: return a full or sparse identity matrix
%-------------------------------------------------------------------------------

function I = identity (A, n)
%IDENTITY return a full or sparse identity matrix.  Not user-callable
if (issparse (A))
I = speye (n) ;
else
I = eye (n) ;
end
end

%-------------------------------------------------------------------------------
% throw an error if inv(A) is being inadvertently computed
%-------------------------------------------------------------------------------

function error_if_inverse (F, kind)
% x = b\F or F/b where F=inverse(A) and b is not a scalar is unsupported.
% It could be done by coercing F into an explicit matrix representation of
% inv(A), via x = b\double(F) or double(A)/b, but this is the same as
% b\inv(A) or inv(A)/b respectively.  That is dangerous, and thus it is
% not done here automatically.
if (F.is_inverse)
if (kind == 1)
s1 = 'B\F' ;
s2 = 'B\double(F)' ;
else
s1 = 'F/B' ;
s2 = 'double(F)/B' ;
end
error ('FACTORIZE:unsupported', ...
['%s where F=inverse(A) requires the explicit computation of the ' ...
'inverse.\nThis is ill-advised, so it is never done automatically.'...
'\nTo force it, use %s instead of %s.\n'], s1, s2, s1) ;
end
end
```