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 % % See also mldivide, lu, chol, ldl, qr, svd. % 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