1classdef factorization 2%FACTORIZATION a generic matrix factorization object 3% Normally, this object is created via the F=factorize(A) function. Users 4% do not need to use this method directly. 5% 6% This is an abstract class that is specialized into 13 different kinds of 7% matrix factorizations: 8% 9% factorization_chol_dense dense Cholesky A = R'*R 10% factorization_lu_dense dense LU A(p,:) = L*U 11% factorization_qr_dense dense QR of A A = Q*R 12% factorization_qrt_dense dense QR of A' A' = Q*R 13% factorization_ldl_dense dense LDL A(p,p) = L*D*L' 14% factorization_cod_dense dense COD A = U*R*V' 15% 16% factorization_chol_sparse sparse Cholesky P'*A*P = L*L' 17% factorization_lu_sparse sparse LU P*(R\A)*Q = L*U 18% factorization_qr_sparse sparse QR of A (A*P)'*(A*P) = R'*R 19% factorization_qrt_sparse sparse QR of A' (P*A)*(P*A)' = R'*R 20% factorization_ldl_sparse sparse LDL P'*A*P = L*D*L' 21% factorization_cod_sparse sparse COD A = U*R*V' 22% 23% factorization_svd SVD A = U*S*V' 24% 25% The abstract class provides the following functions. In the descriptions, 26% F is a factorization. The arguments b, y, and z may be factorizations or 27% matrices. The output x is normally matrix unless it can be represented as a 28% scaled factorization. For example, G=F\2 and G=inverse(F)*2 both return a 29% factorization G. Below, s is always a scalar, and C is always a matrix. 30% 31% These methods return a matrix x, unless one argument is a scalar (in which 32% case they return a scaled factorization object): 33% x = mldivide (F, b) x = A \ b 34% x = mrdivide (b, F) x = b / A 35% x = mtimes (y, z) y * z 36% 37% These methods always return a factorization: 38% F = uplus (F) +F 39% F = uminus (F) -F 40% F = inverse (F) representation of inv(A), without computing it 41% F = ctranspose (F) F' 42% 43% These built-in methods return a scalar: 44% s = isreal (F) 45% s = isempty (F) 46% s = isscalar (F) 47% s = issingle (F) 48% s = isnumeric (F) 49% s = isfloat (F) 50% s = isvector (F) 51% s = issparse (F) 52% s = isfield (F,f) 53% s = isa (F, s) 54% s = condest (F) 55% 56% This method returns the estimated rank from the factorization. 57% s = rankest (F) 58% 59% These methods support access to the contents of a factorization object 60% e = end (F, k, n) 61% [m,n] = size (F, k) 62% S = double (F) 63% C = subsref (F, ij) 64% S = struct (F) 65% disp (F) 66% 67% The factorization_chol_dense object also provides cholupdate, which acts 68% just like the builtin cholupdate. 69% 70% The factorization_svd object provides: 71% 72% c = cond (F,p) the p-norm condition number. p=2 is the default. 73% cond(F,2) takes no time to compute, since it was 74% computed when the SVD factorization was found. 75% a = norm (F,p) the p-norm. see the cond(F,p) discussion above. 76% r = rank (F) returns the rank of A, precomputed by the SVD. 77% Z = null (F) orthonormal basis for the null space of A 78% Q = orth (F) orthonormal basis for the range of A 79% C = pinv (F) the pseudo-inverse, V'*(S\V). 80% [U,S,V] = svd (F) SVD of A or pinv(A), regular, economy, or rank-sized 81% 82% See also mldivide, lu, chol, ldl, qr, svd. 83 84% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com 85 86 properties (SetAccess = protected) 87 % The abstract class holds a QR, LU, Cholesky factorization: 88 A = [ ] ; % a copy of the input matrix 89 Factors = [ ] ; 90 is_inverse = false ;% F represents the factorization of A or inv(A) 91 is_ctrans = false ; % F represents the factorization of A or A' 92 kind = '' ; % a string stating what kind of factorization F is 93 alpha = 1 ; % F represents the factorization of A or alpha*A. 94 A_rank = [ ] ; % rank of A, from SVD, or estimate from COD 95 A_cond = [ ] ; % 2-norm condition number of A, from SVD 96 A_condest = [ ] ; % quick and dirty estimate of the condition number 97 % If F is inverted, alpha doesn't change. For example: 98 % F = alpha*factorize(A) ; % F = alpha*A, in factorized form. 99 % G = inverse(F) ; % G = inv(alpha*A) 100 % H = beta*G % H = beta*inv(alpha*A) 101 % = inv((alpha/beta)*A) 102 % So to update alpha via scaling, beta*F, the new scale factor beta 103 % is multiplied with F.alpha if F.is_inverse is false. Otherwise, 104 % F.alpha is divided by beta to get the new scale factor. 105 end 106 107 methods (Abstract) 108 x = mldivide_subclass (F, b) ; 109 x = mrdivide_subclass (b, F) ; 110 e = error_check (F) ; 111 end 112 113 methods 114 115 %----------------------------------------------------------------------- 116 % mldivide and mrdivide: return a scaled factorization or a matrix 117 %----------------------------------------------------------------------- 118 119 % Let b be a double scalar, F a non-scalar factorization, and g a scalar 120 % factorization. Then these operations return scaled factorization 121 % objects (unless flatten is true, in which case a matrix is returned): 122 % 123 % F\b = inverse (F) * b | F/b = F / b 124 % F\g = inverse (F) * double (g) | F/g = F / double (g) 125 % b\F = F / b | b/F = b * inverse (F) 126 % g\F = F / double (g) | g/F = double (g) * inverse (F) 127 % 128 % Otherwise mldivide & mrdivide always return a matrix as their result. 129 130 function x = mldivide (y, z, flatten) 131 %MLDIVIDE x = y\z where either y or z or both are factorizations. 132 flatten = (nargin > 2 && flatten) ; 133 if (isobject (y) && isscalar (z) && ~flatten) 134 % x = y\z where y is an object and z is scalar (perhaps object). 135 % result is a scaled factorization object x. 136 x = scale_factor (inverse (y), ~(y.is_inverse), double (z)) ; 137 elseif (isscalar (y) && isobject (z) && ~flatten) 138 % x = y\z where y is scalar (perhaps object) and z is an object. 139 % result is a scaled factorization object x. 140 x = scale_factor (z, ~(z.is_inverse), double (y)) ; 141 else 142 % result x will be a matrix. b is coerced to be a matrix. 143 [F, b, first_arg_is_F] = getargs (y, z) ; 144 if (~first_arg_is_F) 145 % x = b\F where F is a factorization. 146 error_if_inverse (F, 1) ; 147 x = b \ F.A ; % use builtin backslash 148 elseif (F.is_ctrans) 149 % x = F\b where F represents (alpha*A)' or inv(alpha*A)' 150 if (F.is_inverse) 151 % x = inv(alpha*A)'\b = (A'*b)*alpha' 152 x = (F.A'*b) * F.alpha' ; 153 else 154 % x = (alpha*A)'\b = (b'/A)' / alpha' 155 x = mrdivide_subclass (b', F)' / F.alpha' ; 156 end 157 else 158 % x = F\b where F represents (alpha*A) or inv(alpha*A) 159 if (F.is_inverse) 160 % x = inv(alpha*A)\b = (A*b)*alpha 161 x = (F.A*b) * F.alpha ; 162 else 163 % x = (alpha*A)\b = (A\b) / alpha 164 x = mldivide_subclass (F, b) / F.alpha ; 165 end 166 end 167 end 168 end 169 170 function x = mrdivide (y, z, flatten) 171 %MRDIVIDE x = y/z where either y or z or both are factorizations. 172 flatten = (nargin > 2 && flatten) ; 173 if (isobject (y) && isscalar (z) && ~flatten) 174 % x = y/z where y is an object and z is scalar (perhaps object). 175 % result is a scaled factorization object x. 176 x = scale_factor (y, ~(y.is_inverse), double (z)) ; 177 elseif (isscalar (y) && isobject (z) && ~flatten) 178 % x = y/z where y is scalar (perhaps object) and z is an object. 179 % result is a scaled factorization object x. 180 x = scale_factor (inverse (z), ~(z.is_inverse), double (y)) ; 181 else 182 % result x will be a matrix. b is coerced to be a matrix. 183 [F, b, first_arg_is_F] = getargs (y, z) ; 184 if (first_arg_is_F) 185 % x = F/b where F is a factorization object 186 error_if_inverse (F, 2) 187 x = F.A / b ; % use builtin slash 188 elseif (F.is_ctrans) 189 % x = b/F where F represents (alpha*A)' or inv(alpha*A)' 190 if (F.is_inverse) 191 % x = b/inv(alpha*A)' = (b*A')*alpha' 192 x = (b*F.A') * F.alpha' ; 193 else 194 % x = b/(alpha*A)' = (A\b')' / alpha' 195 x = mldivide_subclass (F, b')' / F.alpha' ; 196 end 197 else 198 % x = b/F where F represents (alpha*A) or inv(alpha*A) 199 if (F.is_inverse) 200 % x = b/inv(alpha*A) = (b*A)*alpha 201 x = (b*F.A) * F.alpha ; 202 else 203 % x = b/(alpha*A) = (b/A) / alpha 204 x = mrdivide_subclass (b, F) / F.alpha ; 205 end 206 end 207 end 208 end 209 210 %----------------------------------------------------------------------- 211 % mtimes: a simple and clean wrapper for mldivide and mrdivide 212 %----------------------------------------------------------------------- 213 214 function x = mtimes (y, z) 215 %MTIMES x=y*z where y or z is a factorization object (or both). 216 % Since inverse(F) is so cheap, and does the right thing inside 217 % mldivide and mrdivide, this is just a simple wrapper. 218 if (isobject (y)) 219 % A*b becomes inverse(A)\b 220 % inverse(A)*b becomes A\b 221 % A'*b becomes inverse(A)'\b 222 % inverse(A)'*b becomes A'\b 223 x = mldivide (inverse (y), z) ; 224 else 225 % b*A becomes b/inverse(A) 226 % b*inverse(A) becomes b/A 227 % b*A' becomes b/inverse(A)' 228 % b*inverse(A)' becomes b/A' 229 % y is a scalar or matrix, z must be an object 230 x = mrdivide (y, inverse (z)) ; 231 end 232 end 233 234 %----------------------------------------------------------------------- 235 % uplus, uminus, ctranspose, inverse: always return a factorization 236 %----------------------------------------------------------------------- 237 238 function F = uplus (F) 239 %UPLUS +F 240 end 241 242 function F = uminus (F) 243 %UMINUS -F 244 F.alpha = -(F.alpha) ; 245 end 246 247 function F = inverse (F) 248 %INVERSE "inverts" F by flagging it as factorization of inv(A) 249 F.is_inverse = ~(F.is_inverse) ; 250 end 251 252 function F = ctranspose (F) 253 %CTRANSPOSE "transposes" F by flagging it as factorization of A' 254 F.is_ctrans = ~(F.is_ctrans) ; 255 end 256 257 %----------------------------------------------------------------------- 258 % is* methods that return a scalar 259 %----------------------------------------------------------------------- 260 261 function s = isreal (F) 262 %ISREAL for F=factorize(A): same as isreal(A) 263 s = isreal (F.A) ; 264 end 265 266 function s = isempty (F) 267 %ISEMPTY for F=factorize(A): same as isempty(A) 268 s = any (size (F.A) == 0) ; 269 end 270 271 function s = isscalar (F) 272 %ISSCALAR for F=factorize(A): same as isscalar(A) 273 s = isscalar (F.A) ; 274 end 275 276 function s = issingle (F) %#ok 277 %ISSINGLE for F=factorize(A) is always false 278 s = false ; 279 end 280 281 function s = isnumeric (F) %#ok 282 %ISNUMERIC for F=factorize(A) is always true 283 s = true ; 284 end 285 286 function s = isfloat (F) %#ok 287 %ISFLOAT for F=factorize(A) is always true 288 s = true ; 289 end 290 291 function s = isvector (F) 292 %ISVECTOR for F=factorize(A): same as isvector(A) 293 s = isvector (F.A) ; 294 end 295 296 function s = issparse (F) 297 %ISSPARSE for F=factorize(A): same as issparse(A) 298 s = issparse (F.A) ; 299 end 300 301 function s = isfield (F, f) %#ok 302 %ISFIELD isfield(F,f) is true if F.f exists, false otherwise 303 s = (ischar (f) && (strcmp (f, 'A') ... 304 || strcmp (f, 'Factors') || strcmp (f, 'kind') ... 305 || strcmp (f, 'is_inverse') || strcmp (f, 'is_ctrans') ... 306 || strcmp (f, 'alpha') || strcmp (f, 'A_rank') ... 307 || strcmp (f, 'A_cond') || strcmp (f, 'A_condest'))) ; 308 end 309 310 function s = isa (F, s) 311 %ISA for F=factorize(A): 'double', 'numeric', 'float' are true. 312 % For other types, the builtin isa does the right thing. 313 s = strcmp (s, 'double') || strcmp (s, 'numeric') || ... 314 strcmp (s, 'float') || builtin ('isa', F, s) ; 315 end 316 317 %----------------------------------------------------------------------- 318 % condest, rankest 319 %----------------------------------------------------------------------- 320 321 function C = abs (F) 322 %ABS abs(F) returns abs(A) or abs(inverse(A)), as appropriate. The 323 % ONLY reason abs is included here is to support the builtin 324 % normest1 for small matrices (n <= 4). Computing abs(inverse(A)) 325 % explicitly computes the inverse of A, so use with caution. 326 C = abs (double (F)) ; 327 end 328 329 function s = condest (F) 330 %CONDEST 1-norm condition number for square matrices. 331 % Does not require another factorization of A, so it's very fast. 332 % Does NOT explicitly compute the inverse of A. Instead, if F 333 % represents an inverse, F*x inside normest1 does the right thing, 334 % and does A\b using the factorization F. 335 A = F.A ; %#ok 336 [m, n] = size (A) ; %#ok 337 if (m ~= n) 338 error ('MATLAB:condest:NonSquareMatrix', ... 339 'Matrix must be square.') ; 340 end 341 if (n == 0) 342 s = 0 ; 343 elseif (F.is_inverse) 344 % F already represents the factorization of the inverse of A 345 s = F.alpha * norm (A,1) * normest1 (F) ; %#ok 346 else 347 % Note that the inverse is NOT explicitly computed. 348 s = F.alpha * norm (A,1) * normest1 (inverse (F)) ; %#ok 349 end 350 end 351 352 function r = rankest (F) 353 %RANKEST returns the estimated rank of A. 354 % It is a very rough estimate if Cholesky, LU, QR, or LDL succeeded 355 % (in which A is assumed to have full rank). COD finds a more 356 % accurate estimate, and SVD finds the exact rank. 357 r = F.A_rank ; 358 end 359 360 %----------------------------------------------------------------------- 361 % end, size 362 %----------------------------------------------------------------------- 363 364 function e = end (F, k, n) 365 %END returns index of last item for use in subsref 366 if (n == 1) 367 e = numel (F.A) ; % # of elements, for linear indexing 368 else 369 e = size (F, k) ; % # of rows or columns in A or pinv(A) 370 end 371 end 372 373 function [m, n] = size (F, k) 374 %SIZE returns the size of the matrix F.A in the factorization F 375 if (F.is_inverse ~= F.is_ctrans) 376 % swap the dimensions to match pinv(A) 377 if (nargout > 1) 378 [n, m] = size (F.A) ; 379 else 380 m = size (F.A) ; 381 m = m ([2 1]) ; 382 end 383 else 384 if (nargout > 1) 385 [m, n] = size (F.A) ; 386 else 387 m = size (F.A) ; 388 end 389 end 390 if (nargin > 1) 391 m = m (k) ; 392 end 393 end 394 395 %----------------------------------------------------------------------- 396 % double: a wrapper for subsref 397 %----------------------------------------------------------------------- 398 399 function S = double (F) 400 %DOUBLE returns the factorization as a matrix, A or inv(A) 401 ij.type = '()' ; 402 ij.subs = cell (1,0) ; 403 S = subsref (F, ij) ; % let factorize.subsref do all the work 404 end 405 406 %----------------------------------------------------------------------- 407 % subsref: returns a matrix 408 %----------------------------------------------------------------------- 409 410 function C = subsref (F, ij) 411 %SUBSREF A(i,j) or (i,j)th entry of inv(A) if F is inverted. 412 % Otherwise, explicit entries in the inverse are computed. 413 % This method also extracts the contents of F with F.whatever. 414 switch (ij (1).type) 415 case '.' 416 % F.A usage: extract one of the matrices from F 417 switch ij (1).subs 418 case 'A' 419 C = F.A ; 420 case 'Factors' 421 C = F.Factors ; 422 case 'is_inverse' 423 C = F.is_inverse ; 424 case 'is_ctrans' 425 C = F.is_ctrans ; 426 case 'kind' 427 C = F.kind ; 428 case 'alpha' 429 C = F.alpha ; 430 case 'A_cond' 431 C = F.A_cond ; 432 case 'A_condest' 433 C = F.A_condest ; 434 case 'A_rank' 435 C = F.A_rank ; 436 otherwise 437 error ('MATLAB:nonExistentField', ... 438 'Reference to non-existent field ''%s''.', ... 439 ij (1).subs) ; 440 end 441 % F.X(2,3) usage, return X(2,3), for component F.X 442 if (length (ij) > 1 && ~isempty (ij (2).subs)) 443 C = subsref (C, ij (2)) ; 444 end 445 case '()' 446 C = subsref_paren (F, ij) ; 447 case '{}' 448 error ('MATLAB:cellRefFromNonCell', ... 449 'Cell contents reference from a non-cell array object.') ; 450 end 451 end 452 453 %----------------------------------------------------------------------- 454 % struct: extracts all contents of a factorization object 455 %----------------------------------------------------------------------- 456 457 function S = struct (F) 458 %STRUCT convert factorization F into a struct. 459 % S cannot be used for subsequent object methods here. 460 S.A = F.A ; 461 S.Factors = F.Factors ; 462 S.is_inverse = F.is_inverse ; 463 S.is_ctrans = F.is_ctrans ; 464 S.alpha = F.alpha ; 465 S.A_rank = F.A_rank ; 466 S.A_cond = F.A_cond ; 467 S.kind = F.kind ; 468 end 469 470 %----------------------------------------------------------------------- 471 % disp: displays the contents of F 472 %----------------------------------------------------------------------- 473 474 function disp (F) 475 %DISP displays a factorization object 476 fprintf (' class: %s\n', class (F)) ; 477 fprintf (' %s\n', F.kind) ; 478 fprintf (' A: [%dx%d double]\n', size (F.A)) ; 479 fprintf (' Factors:\n') ; disp (F.Factors) ; 480 fprintf (' is_inverse: %d\n', F.is_inverse) ; 481 fprintf (' is_ctrans: %d\n', F.is_ctrans) ; 482 fprintf (' alpha: %g', F.alpha) ; 483 if (~isreal (F.alpha)) 484 fprintf (' + (%g)i', imag (F.alpha)) ; 485 end 486 fprintf ('\n') ; 487 if (~isempty (F.A_rank)) 488 fprintf (' A_rank: %d\n', F.A_rank) ; 489 end 490 if (~isempty (F.A_condest)) 491 fprintf (' A_condest: %d\n', F.A_condest) ; 492 end 493 if (~isempty (F.A_cond)) 494 fprintf (' A_cond: %d\n', F.A_cond) ; 495 end 496 end 497 end 498 499 %--------------------------------------------------------------------------- 500 % methods that are not user-callable 501 %--------------------------------------------------------------------------- 502 503 methods (Access = protected) 504 505 function [F, b, first_arg_is_F] = getargs (y, z) 506 first_arg_is_F = isobject (y) ; 507 if (first_arg_is_F) 508 F = y ; % first argument is a factorization object 509 b = double (z) ; % 2nd one coerced to be a matrix 510 else 511 b = y ; % first argument is not an object 512 F = z ; % second one must be an object 513 end 514 end 515 516 function F = scale_factor (F, use_beta_inverse, beta) 517 %SCALE_FACTOR scales a factorization 518 if (use_beta_inverse) 519 % F = inv(alpha*A), so F*beta = inv((alpha/beta)*A) 520 if (F.is_ctrans) 521 F.alpha = F.alpha / beta' ; 522 else 523 F.alpha = F.alpha / beta ; 524 end 525 else 526 % F = alpha*A, so F*beta = (alpha*beta)*A 527 if (F.is_ctrans) 528 F.alpha = F.alpha * beta' ; 529 else 530 F.alpha = F.alpha * beta ; 531 end 532 end 533 end 534 end 535end 536 537 538%------------------------------------------------------------------------------- 539% subsref_paren: support function for subsref 540%------------------------------------------------------------------------------- 541 542function C = subsref_paren (F, ij) 543%SUBSREF_PAREN C = subsref_paren(F,ij) implements C=F(i,j) and C=F(i) 544 545 % F(2,3) usage, return A(2,3) or the (2,3) of inv(A). 546 assert (length (ij) == 1, 'Improper index matrix reference.') ; 547 A = F.A ; 548 is_ctrans = F.is_ctrans ; 549 if (is_ctrans && length (ij.subs) > 1) % swap i and j 550 ij.subs = ij.subs ([2 1]) ; 551 end 552 553 if (F.is_inverse) 554 555 % requesting explicit entries of the inverse 556 557 if (length (ij.subs) == 1) 558 % for linear indexing of the inverse (C=F(i)), first 559 % convert to double and then use builtin subsref 560 C = subsref (double (F), ij) ; 561 else 562 % standard indexing, C = F(i,j) 563 if (is_ctrans) 564 [n, m] = size (A) ; 565 else 566 [m, n] = size (A) ; 567 end 568 if (length (ij.subs) == 2) 569 ilen = length (ij.subs {1}) ; 570 if (strcmp (ij.subs {1}, ':')) 571 ilen = n ; 572 end 573 jlen = length (ij.subs {2}) ; 574 if (strcmp (ij.subs {2}, ':')) 575 jlen = m ; 576 end 577 j = ij ; 578 j.subs {1} = ':' ; 579 i = ij ; 580 i.subs {2} = ':' ; 581 if (jlen <= ilen) 582 % compute X=S(:,j) of S=inv(A) and return X(i,:) 583 C = subsref (mldivide (... 584 inverse (F), ... 585 subsref (identity (A, m), j), 1), i) ; 586 else 587 % compute X=S(i,:) of S=inv(A) and return X(:,j) 588 C = subsref (mrdivide (... 589 subsref (identity (A, n), i), ... 590 inverse (F), 1), j) ; 591 end 592 else 593 % the entire inverse has been explicitly computed 594 C = mldivide (inverse (F), identity (A, m), 1) ; 595 end 596 end 597 598 else 599 600 % F is not inverted, so just return A(i,j) 601 if (isempty (ij (1).subs)) 602 C = A ; 603 else 604 C = subsref (A, ij) ; 605 end 606 C = C * F.alpha ; 607 if (is_ctrans) 608 C = C' ; 609 end 610 end 611end 612 613 614%------------------------------------------------------------------------------- 615% identity: return a full or sparse identity matrix 616%------------------------------------------------------------------------------- 617 618function I = identity (A, n) 619%IDENTITY return a full or sparse identity matrix. Not user-callable 620 if (issparse (A)) 621 I = speye (n) ; 622 else 623 I = eye (n) ; 624 end 625end 626 627%------------------------------------------------------------------------------- 628% throw an error if inv(A) is being inadvertently computed 629%------------------------------------------------------------------------------- 630 631function error_if_inverse (F, kind) 632 % x = b\F or F/b where F=inverse(A) and b is not a scalar is unsupported. 633 % It could be done by coercing F into an explicit matrix representation of 634 % inv(A), via x = b\double(F) or double(A)/b, but this is the same as 635 % b\inv(A) or inv(A)/b respectively. That is dangerous, and thus it is 636 % not done here automatically. 637 if (F.is_inverse) 638 if (kind == 1) 639 s1 = 'B\F' ; 640 s2 = 'B\double(F)' ; 641 else 642 s1 = 'F/B' ; 643 s2 = 'double(F)/B' ; 644 end 645 error ('FACTORIZE:unsupported', ... 646 ['%s where F=inverse(A) requires the explicit computation of the ' ... 647 'inverse.\nThis is ill-advised, so it is never done automatically.'... 648 '\nTo force it, use %s instead of %s.\n'], s1, s2, s1) ; 649 end 650end 651