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 p*A = 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 q'*A*q = L*D*L'
% factorization_cod_dense dense COD A = U*R*V'
%
% factorization_chol_sparse sparse Cholesky q'*A*q = L*L'
% factorization_lu_sparse sparse LU p*A*q = L*U
% factorization_qr_sparse sparse QR of A (A*q)'*(A*q) = R'*R
% factorization_qrt_sparse sparse QR of A' (p*A)*(p*A)' = R'*R
% factorization_ldl_sparse sparse LDL q'*A*q = 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. s is always a scalar. 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 methods return matrices / scalars
% [i,j,x] = find (F,...)
% s = nnz (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
%
% See also mldivide, lu, chol, ldl, qr, svd.
% Copyright 2011, Timothy A. Davis, University of Florida.
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) ;
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
