octave_packages/linear-algebra-2.2.0/doc-cache

   1 # Created by Octave 3.6.1, Mon May 21 19:48:13 2012 UTC <root@brouzouf>
   2 # name: cache
   3 # type: cell
   4 # rows: 3
   5 # columns: 16
   6 # name: <cell-element>
   7 # type: sq_string
   8 # elements: 1
   9 # length: 8
  10 cartprod
  11
  12
  13 # name: <cell-element>
  14 # type: sq_string
  15 # elements: 1
  16 # length: 770
  17  -- Function File:  cartprod (VARARGIN)
  18      Computes the cartesian product of given column vectors ( row
  19      vectors ).  The vector elements are assumend to be numbers.
  20
  21      Alternatively the vectors can be specified by as a matrix, by its
  22      columns.
  23
  24      To calculate the cartesian product of vectors, P = A x B x C x D
  25      ... . Requires A, B, C, D be column vectors.  The algorithm is
  26      iteratively calcualte the products,  ( ( (A x B ) x C ) x D ) x
  27      etc.
  28
  29             cartprod(1:2,3:4,0:1)
  30             ans =   1   3   0
  31                     2   3   0
  32                     1   4   0
  33                     2   4   0
  34                     1   3   1
  35                     2   3   1
  36                     1   4   1
  37                     2   4   1
  38
  39    See also: kron
  40
  41
  42
  43
  44 # name: <cell-element>
  45 # type: sq_string
  46 # elements: 1
  47 # length: 71
  48 Computes the cartesian product of given column vectors ( row vectors ).
  49
  50
  51
  52 # name: <cell-element>
  53 # type: sq_string
  54 # elements: 1
  55 # length: 13
  56 circulant_eig
  57
  58
  59 # name: <cell-element>
  60 # type: sq_string
  61 # elements: 1
  62 # length: 718
  63  -- Function File: LAMBDA = circulant_eig (V)
  64  -- Function File: [VS, LAMBDA] = circulant_eig (V)
  65      Fast, compact calculation of eigenvalues and eigenvectors of a
  66      circulant matrix
  67      Given an N*1 vector V, return the eigenvalues LAMBDA and
  68      optionally eigenvectors VS of the N*N circulant matrix C that has
  69      V as its first column
  70
  71      Theoretically same as `eig(make_circulant_matrix(v))', but many
  72      fewer computations; does not form C explicitly
  73
  74      Reference: Robert M. Gray, Toeplitz and Circulant Matrices: A
  75      Review, Now Publishers, http://ee.stanford.edu/~gray/toeplitz.pdf,
  76      Chapter 3
  77
  78      See also: circulant_make_matrix, circulant_matrix_vector_product,
  79      circulant_inv
  80
  81
  82
  83
  84
  85 # name: <cell-element>
  86 # type: sq_string
  87 # elements: 1
  88 # length: 80
  89 Fast, compact calculation of eigenvalues and eigenvectors of a
  90 circulant matrix
  91
  92
  93
  94
  95 # name: <cell-element>
  96 # type: sq_string
  97 # elements: 1
  98 # length: 13
  99 circulant_inv
 100
 101
 102 # name: <cell-element>
 103 # type: sq_string
 104 # elements: 1
 105 # length: 875
 106  -- Function File: C = circulant_inv (V)
 107      Fast, compact calculation of inverse of a circulant matrix
 108      Given an N*1 vector V, return the inverse C of the N*N circulant
 109      matrix C that has V as its first column The returned C is the
 110      first column of the inverse, which is also circulant - to get the
 111      full matrix, use `circulant_make_matrix(c)'
 112
 113      Theoretically same as `inv(make_circulant_matrix(v))(:, 1)', but
 114      requires many fewer computations and does not form matrices
 115      explicitly
 116
 117      Roundoff may induce a small imaginary component in C even if V is
 118      real - use `real(c)' to remedy this
 119
 120      Reference: Robert M. Gray, Toeplitz and Circulant Matrices: A
 121      Review, Now Publishers, http://ee.stanford.edu/~gray/toeplitz.pdf,
 122      Chapter 3
 123
 124      See also: circulant_make_matrix, circulant_matrix_vector_product,
 125      circulant_eig
 126
 127
 128
 129
 130
 131 # name: <cell-element>
 132 # type: sq_string
 133 # elements: 1
 134 # length: 80
 135 Fast, compact calculation of inverse of a circulant matrix
 136 Given an N*1 vector V
 137
 138
 139
 140 # name: <cell-element>
 141 # type: sq_string
 142 # elements: 1
 143 # length: 21
 144 circulant_make_matrix
 145
 146
 147 # name: <cell-element>
 148 # type: sq_string
 149 # elements: 1
 150 # length: 579
 151  -- Function File: C = circulant_make_matrix (V)
 152      Produce a full circulant matrix given the first column
 153      Given an N*1 vector V, returns the N*N circulant matrix C where V
 154      is the left column and all other columns are downshifted versions
 155      of V
 156
 157      Note: If the first row R of a circulant matrix is given, the first
 158      column V can be obtained as `v = r([1 end:-1:2])'
 159
 160      Reference: Gene H. Golub and Charles F. Van Loan, Matrix
 161      Computations, 3rd Ed., Section 4.7.7
 162
 163      See also: circulant_matrix_vector_product, circulant_eig,
 164      circulant_inv
 165
 166
 167
 168
 169
 170 # name: <cell-element>
 171 # type: sq_string
 172 # elements: 1
 173 # length: 80
 174 Produce a full circulant matrix given the first column
 175 Given an N*1 vector V, re
 176
 177
 178
 179 # name: <cell-element>
 180 # type: sq_string
 181 # elements: 1
 182 # length: 31
 183 circulant_matrix_vector_product
 184
 185
 186 # name: <cell-element>
 187 # type: sq_string
 188 # elements: 1
 189 # length: 751
 190  -- Function File: Y = circulant_matrix_vector_product (V, X)
 191      Fast, compact calculation of the product of a circulant matrix
 192      with a vector
 193      Given N*1 vectors V and X, return the matrix-vector product Y =
 194      CX, where C is the N*N circulant matrix that has V as its first
 195      column
 196
 197      Theoretically the same as `make_circulant_matrix(x) * v', but does
 198      not form C explicitly; uses the discrete Fourier transform
 199
 200      Because of roundoff, the returned Y may have a small imaginary
 201      component even if V and X are real (use `real(y)' to remedy this)
 202
 203      Reference: Gene H. Golub and Charles F. Van Loan, Matrix
 204      Computations, 3rd Ed., Section 4.7.7
 205
 206      See also: circulant_make_matrix, circulant_eig, circulant_inv
 207
 208
 209
 210
 211
 212 # name: <cell-element>
 213 # type: sq_string
 214 # elements: 1
 215 # length: 80
 216 Fast, compact calculation of the product of a circulant matrix with a
 217 vector
 218 Giv
 219
 220
 221
 222 # name: <cell-element>
 223 # type: sq_string
 224 # elements: 1
 225 # length: 3
 226 cod
 227
 228
 229 # name: <cell-element>
 230 # type: sq_string
 231 # elements: 1
 232 # length: 887
 233  -- Function File: [Q, R, Z] = cod (A)
 234  -- Function File: [Q, R, Z, P] = cod (A)
 235  -- Function File: [...] = cod (A, '0')
 236      Computes the complete orthogonal decomposition (COD) of the matrix
 237      A:
 238             A = Q*R*Z'
 239      Let A be an M-by-N matrix, and let `K = min(M, N)'.  Then Q is
 240      M-by-M orthogonal, Z is N-by-N orthogonal, and R is M-by-N such
 241      that `R(:,1:K)' is upper trapezoidal and `R(:,K+1:N)' is zero.
 242      The additional P output argument specifies that pivoting should be
 243      used in the first step (QR decomposition). In this case,
 244             A*P = Q*R*Z'
 245      If a second argument of '0' is given, an economy-sized
 246      factorization is returned so that R is K-by-K.
 247
 248      _NOTE_: This is currently implemented by double QR factorization
 249      plus some tricky manipulations, and is not as efficient as using
 250      xRZTZF from LAPACK.
 251
 252      See also: qr
 253
 254
 255
 256
 257
 258 # name: <cell-element>
 259 # type: sq_string
 260 # elements: 1
 261 # length: 80
 262 Computes the complete orthogonal decomposition (COD) of the matrix A:
 263        A =
 264
 265
 266
 267 # name: <cell-element>
 268 # type: sq_string
 269 # elements: 1
 270 # length: 7
 271 condeig
 272
 273
 274 # name: <cell-element>
 275 # type: sq_string
 276 # elements: 1
 277 # length: 784
 278  -- Function File: C = condeig (A)
 279  -- Function File: [V, LAMBDA, C] = condeig (A)
 280      Compute condition numbers of the eigenvalues of a matrix. The
 281      condition numbers are the reciprocals of the cosines of the angles
 282      between the left and right eigenvectors.
 283
 284 Arguments
 285 ---------
 286
 287         * A must be a square numeric matrix.
 288
 289 Return values
 290 -------------
 291
 292         * C is a vector of condition numbers of the eigenvalue of A.
 293
 294         * V is the matrix of right eigenvectors of A. The result is the
 295           same as for `[v, lambda] = eig (a)'.
 296
 297         * LAMBDA is the diagonal matrix of eigenvalues of A. The result
 298           is the same as for `[v, lambda] = eig (a)'.
 299
 300 Example
 301 -------
 302
 303           a = [1, 2; 3, 4];
 304           c = condeig (a)
 305           => [1.0150; 1.0150]
 306
 307
 308
 309
 310 # name: <cell-element>
 311 # type: sq_string
 312 # elements: 1
 313 # length: 57
 314 Compute condition numbers of the eigenvalues of a matrix.
 315
 316
 317
 318 # name: <cell-element>
 319 # type: sq_string
 320 # elements: 1
 321 # length: 4
 322 funm
 323
 324
 325 # name: <cell-element>
 326 # type: sq_string
 327 # elements: 1
 328 # length: 1117
 329  -- Function File: B = funm (A, F)
 330      Compute matrix equivalent of function F; F can be a function name
 331      or a function handle.
 332
 333      For trigonometric and hyperbolic functions, `thfm' is automatically
 334      invoked as that is based on `expm' and diagonalization is avoided.
 335      For other functions diagonalization is invoked, which implies that
 336      -depending on the properties of input matrix A- the results can be
 337      very inaccurate _without any warning_. For easy diagonizable and
 338      stable matrices results of funm will be sufficiently accurate.
 339
 340      Note that you should not use funm for 'sqrt', 'log' or 'exp';
 341      instead use sqrtm, logm and expm as these are more robust.
 342
 343      Examples:
 344
 345             B = funm (A, sin);
 346             (Compute matrix equivalent of sin() )
 347
 348             function bk1 = besselk1 (x)
 349                bk1 = besselk(x, 1);
 350             endfunction
 351             B = funm (A, besselk1);
 352             (Compute matrix equivalent of bessel function K1(); a helper function
 353              is needed here to convey extra args for besselk() )
 354
 355      See also: thfm, expm, logm, sqrtm
 356
 357
 358
 359
 360
 361 # name: <cell-element>
 362 # type: sq_string
 363 # elements: 1
 364 # length: 80
 365 Compute matrix equivalent of function F; F can be a function name or a
 366 function
 367
 368
 369
 370 # name: <cell-element>
 371 # type: sq_string
 372 # elements: 1
 373 # length: 6
 374 lobpcg
 375
 376
 377 # name: <cell-element>
 378 # type: sq_string
 379 # elements: 1
 380 # length: 9717
 381  -- Function File: [BLOCKVECTORX, LAMBDA] = lobpcg (BLOCKVECTORX,
 382           OPERATORA)
 383  -- Function File: [BLOCKVECTORX, LAMBDA, FAILUREFLAG] = lobpcg
 384           (BLOCKVECTORX, OPERATORA)
 385  -- Function File: [BLOCKVECTORX, LAMBDA, FAILUREFLAG, LAMBDAHISTORY,
 386 RESIDUALNORMSHISTORY] = lobpcg (BLOCKVECTORX, OPERATORA, OPERATORB,
 387           OPERATORT, BLOCKVECTORY, RESIDUALTOLERANCE, MAXITERATIONS,
 388           VERBOSITYLEVEL)
 389      Solves Hermitian partial eigenproblems using preconditioning.
 390
 391      The first form outputs the array of algebraic smallest eigenvalues
 392      LAMBDA and corresponding matrix of orthonormalized eigenvectors
 393      BLOCKVECTORX of the Hermitian (full or sparse) operator OPERATORA
 394      using input matrix BLOCKVECTORX as an initial guess, without
 395      preconditioning, somewhat similar to:
 396
 397           # for real symmetric operator operatorA
 398           opts.issym  = 1; opts.isreal = 1; K = size (blockVectorX, 2);
 399           [blockVectorX, lambda] = eigs (operatorA, K, 'SR', opts);
 400
 401           # for Hermitian operator operatorA
 402           K = size (blockVectorX, 2);
 403           [blockVectorX, lambda] = eigs (operatorA, K, 'SR');
 404
 405      The second form returns a convergence flag. If FAILUREFLAG is 0
 406      then all the eigenvalues converged; otherwise not all converged.
 407
 408      The third form computes smallest eigenvalues LAMBDA and
 409      corresponding eigenvectors BLOCKVECTORX of the generalized
 410      eigenproblem Ax=lambda Bx, where Hermitian operators OPERATORA and
 411      OPERATORB are given as functions, as well as a preconditioner,
 412      OPERATORT. The operators OPERATORB and OPERATORT must be in
 413      addition _positive definite_. To compute the largest eigenpairs of
 414      OPERATORA, simply apply the code to OPERATORA multiplied by -1.
 415      The code does not involve _any_ matrix factorizations of OPERATORA
 416      and OPERATORB, thus, e.g., it preserves the sparsity and the
 417      structure of OPERATORA and OPERATORB.
 418
 419      RESIDUALTOLERANCE and MAXITERATIONS control tolerance and max
 420      number of steps, and VERBOSITYLEVEL = 0, 1, or 2 controls the
 421      amount of printed info. LAMBDAHISTORY is a matrix with all
 422      iterative lambdas, and RESIDUALNORMSHISTORY are matrices of the
 423      history of 2-norms of residuals
 424
 425      Required input:
 426         * BLOCKVECTORX (class numeric) - initial approximation to
 427           eigenvectors, full or sparse matrix n-by-blockSize.
 428           BLOCKVECTORX must be full rank.
 429
 430         * OPERATORA (class numeric, char, or function_handle) - the
 431           main operator of the eigenproblem, can be a matrix, a
 432           function name, or handle
 433
 434      Optional function input:
 435         * OPERATORB (class numeric, char, or function_handle) - the
 436           second operator, if solving a generalized eigenproblem, can
 437           be a matrix, a function name, or handle; by default if empty,
 438           `operatorB = I'.
 439
 440         * OPERATORT (class char or function_handle) - the
 441           preconditioner, by default `operatorT(blockVectorX) =
 442           blockVectorX'.
 443
 444      Optional constraints input:
 445         * BLOCKVECTORY (class numeric) - a full or sparse n-by-sizeY
 446           matrix of constraints, where sizeY < n. BLOCKVECTORY must be
 447           full rank. The iterations will be performed in the
 448           (operatorB-) orthogonal complement of the column-space of
 449           BLOCKVECTORY.
 450
 451      Optional scalar input parameters:
 452         * RESIDUALTOLERANCE (class numeric) - tolerance, by default,
 453           `residualTolerance = n * sqrt (eps)'
 454
 455         * MAXITERATIONS - max number of iterations, by default,
 456           `maxIterations = min (n, 20)'
 457
 458         * VERBOSITYLEVEL - either 0 (no info), 1, or 2 (with pictures);
 459           by default, `verbosityLevel = 0'.
 460
 461      Required output:
 462         * BLOCKVECTORX and LAMBDA (class numeric) both are computed
 463           blockSize eigenpairs, where `blockSize = size (blockVectorX,
 464           2)' for the initial guess BLOCKVECTORX if it is full rank.
 465
 466      Optional output:
 467         * FAILUREFLAG (class integer) as described above.
 468
 469         * LAMBDAHISTORY (class numeric) as described above.
 470
 471         * RESIDUALNORMSHISTORY (class numeric) as described above.
 472
 473      Functions `operatorA(blockVectorX)', `operatorB(blockVectorX)' and
 474      `operatorT(blockVectorX)' must support BLOCKVECTORX being a
 475      matrix, not just a column vector.
 476
 477      Every iteration involves one application of OPERATORA and
 478      OPERATORB, and one of OPERATORT.
 479
 480      Main memory requirements: 6 (9 if `isempty(operatorB)=0') matrices
 481      of the same size as BLOCKVECTORX, 2 matrices of the same size as
 482      BLOCKVECTORY (if present), and two square matrices of the size
 483      3*blockSize.
 484
 485      In all examples below, we use the Laplacian operator in a 20x20
 486      square with the mesh size 1 which can be generated in MATLAB by
 487      running:
 488           A = delsq (numgrid ('S', 21));
 489           n = size (A, 1);
 490
 491      or in MATLAB and Octave by:
 492           [~,~,A] = laplacian ([19, 19]);
 493           n = size (A, 1);
 494
 495      Note that `laplacian' is a function of the specfun octave-forge
 496      package.
 497
 498      The following Example:
 499           [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, 1e-5, 50, 2);
 500
 501      attempts to compute 8 first eigenpairs without preconditioning,
 502      but not all eigenpairs converge after 50 steps, so failureFlag=1.
 503
 504      The next Example:
 505           blockVectorY = [];
 506           lambda_all = [];
 507           for j = 1:4
 508             [blockVectorX, lambda] = lobpcg (randn (n, 2), A, blockVectorY, 1e-5, 200, 2);
 509             blockVectorY           = [blockVectorY, blockVectorX];
 510             lambda_all             = [lambda_all' lambda']';
 511             pause;
 512           end
 513
 514      attemps to compute the same 8 eigenpairs by calling the code 4
 515      times with blockSize=2 using orthogonalization to the previously
 516      founded eigenvectors.
 517
 518      The following Example:
 519           R       = ichol (A, struct('michol', 'on'));
 520           precfun = @(x)R\(R'\x);
 521           [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, [], @(x)precfun(x), 1e-5, 60, 2);
 522
 523      computes the same eigenpairs in less then 25 steps, so that
 524      failureFlag=0 using the preconditioner function `precfun', defined
 525      inline. If `precfun' is defined as an octave function in a file,
 526      the function handle `@(x)precfun(x)' can be equivalently replaced
 527      by the function name `precfun'. Running:
 528
 529           [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, speye (n), @(x)precfun(x), 1e-5, 50, 2);
 530
 531      produces similar answers, but is somewhat slower and needs more
 532      memory as technically a generalized eigenproblem with B=I is
 533      solved here.
 534
 535      The following example for a mostly diagonally dominant sparse
 536      matrix A demonstrates different types of preconditioning, compared
 537      to the standard use of the main diagonal of A:
 538
 539           clear all; close all;
 540           n       = 1000;
 541           M       = spdiags ([1:n]', 0, n, n);
 542           precfun = @(x)M\x;
 543           A       = M + sprandsym (n, .1);
 544           Xini    = randn (n, 5);
 545           maxiter = 15;
 546           tol     = 1e-5;
 547           [~,~,~,~,rnp] = lobpcg (Xini, A, tol, maxiter, 1);
 548           [~,~,~,~,r]   = lobpcg (Xini, A, [], @(x)precfun(x), tol, maxiter, 1);
 549           subplot (2,2,1), semilogy (r'); hold on;
 550           semilogy (rnp', ':>');
 551           title ('No preconditioning (top)'); axis tight;
 552           M(1,2)  = 2;
 553           precfun = @(x)M\x; % M is no longer symmetric
 554           [~,~,~,~,rns] = lobpcg (Xini, A, [], @(x)precfun(x), tol, maxiter, 1);
 555           subplot (2,2,2), semilogy (r'); hold on;
 556           semilogy (rns', '--s');
 557           title ('Nonsymmetric preconditioning (square)'); axis tight;
 558           M(1,2)  = 0;
 559           precfun = @(x)M\(x+10*sin(x)); % nonlinear preconditioning
 560           [~,~,~,~,rnl] = lobpcg (Xini, A, [], @(x)precfun(x), tol, maxiter, 1);
 561           subplot (2,2,3),  semilogy (r'); hold on;
 562           semilogy (rnl', '-.*');
 563           title ('Nonlinear preconditioning (star)'); axis tight;
 564           M       = abs (M - 3.5 * speye (n, n));
 565           precfun = @(x)M\x;
 566           [~,~,~,~,rs] = lobpcg (Xini, A, [], @(x)precfun(x), tol, maxiter, 1);
 567           subplot (2,2,4),  semilogy (r'); hold on;
 568           semilogy (rs', '-d');
 569           title ('Selective preconditioning (diamond)'); axis tight;
 570
 571 References
 572 ==========
 573
 574      This main function `lobpcg' is a version of the preconditioned
 575 conjugate gradient method (Algorithm 5.1) described in A. V. Knyazev,
 576 Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block
 577 Preconditioned Conjugate Gradient Method, SIAM Journal on Scientific
 578 Computing 23 (2001), no. 2, pp. 517-541.
 579 `http://dx.doi.org/10.1137/S1064827500366124'
 580
 581 Known bugs/features
 582 ===================
 583
 584         * an excessively small requested tolerance may result in often
 585           restarts and instability. The code is not written to produce
 586           an eps-level accuracy! Use common sense.
 587
 588         * the code may be very sensitive to the number of eigenpairs
 589           computed, if there is a cluster of eigenvalues not completely
 590           included, cf.
 591                operatorA = diag ([1 1.99 2:99]);
 592                [blockVectorX, lambda] = lobpcg (randn (100, 1),operatorA, 1e-10, 80, 2);
 593                [blockVectorX, lambda] = lobpcg (randn (100, 2),operatorA, 1e-10, 80, 2);
 594                [blockVectorX, lambda] = lobpcg (randn (100, 3),operatorA, 1e-10, 80, 2);
 595
 596 Distribution
 597 ============
 598
 599      The main distribution site: `http://math.ucdenver.edu/~aknyazev/'
 600
 601      A C-version of this code is a part of the
 602 `http://code.google.com/p/blopex/' package and is directly available,
 603 e.g., in PETSc and HYPRE.
 604
 605
 606
 607
 608 # name: <cell-element>
 609 # type: sq_string
 610 # elements: 1
 611 # length: 61
 612 Solves Hermitian partial eigenproblems using preconditioning.
 613
 614
 615
 616 # name: <cell-element>
 617 # type: sq_string
 618 # elements: 1
 619 # length: 7
 620 ndcovlt
 621
 622
 623 # name: <cell-element>
 624 # type: sq_string
 625 # elements: 1
 626 # length: 771
 627  -- Function File: Y = ndcovlt (X, T1, T2, ...)
 628      Computes an n-dimensional covariant linear transform of an n-d
 629      tensor, given a transformation matrix for each dimension. The
 630      number of columns of each transformation matrix must match the
 631      corresponding extent of X, and the number of rows determines the
 632      corresponding extent of Y. For example:
 633
 634             size (X, 2) == columns (T2)
 635             size (Y, 2) == rows (T2)
 636
 637      The element `Y(i1, i2, ...)' is defined as a sum of
 638
 639             X(j1, j2, ...) * T1(i1, j1) * T2(i2, j2) * ...
 640
 641      over all j1, j2, .... For two dimensions, this reduces to
 642             Y = T1 * X * T2.'
 643
 644      [] passed as a transformation matrix is converted to identity
 645      matrix for the corresponding dimension.
 646
 647
 648
 649
 650
 651 # name: <cell-element>
 652 # type: sq_string
 653 # elements: 1
 654 # length: 80
 655 Computes an n-dimensional covariant linear transform of an n-d tensor,
 656 given a t
 657
 658
 659
 660 # name: <cell-element>
 661 # type: sq_string
 662 # elements: 1
 663 # length: 8
 664 nmf_bpas
 665
 666
 667 # name: <cell-element>
 668 # type: sq_string
 669 # elements: 1
 670 # length: 3768
 671  -- Function File: [W, H, ITER, HIS] =  nmf_bpas (A, K)
 672      Nonnegative Matrix Factorization by Alternating Nonnegativity
 673      Constrained Least Squares using Block Principal Pivoting/Active
 674      Set method.
 675
 676      This function solves one the following problems: given A and K,
 677      find W and H such that     (1) minimize 1/2 * || A-WH ||_F^2
 678      (2) minimize 1/2 * ( || A-WH ||_F^2 + alpha * || W ||_F^2 + beta *
 679      || H ||_F^2 )     (3) minimize 1/2 * ( || A-WH ||_F^2 + alpha * ||
 680      W ||_F^2 + beta * (sum_(i=1)^n || H(:,i) ||_1^2 ) )     where W>=0
 681      and H>=0 elementwise.  The input arguments are A : Input data
 682      matrix (m x n) and K : Target low-rank.
 683
 684      *Optional Inputs*
 685     `Type : Default is 'regularized', which is recommended for quick application testing unless 'sparse' or 'plain' is explicitly needed. If sparsity is needed for 'W' factor, then apply this function for the transpose of 'A' with formulation (3). Then, exchange 'W' and 'H' and obtain the transpose of them. Imposing sparsity for both factors is not recommended and thus not included in this software.'
 686
 687          'plain' to use formulation (1)
 688
 689          'regularized' to use formulation (2)
 690
 691          'sparse' to use formulation (3)
 692
 693     `NNLSSolver : Default is 'bp', which is in general faster.'
 694                item 'bp' to use the algorithm in [1] item 'as' to use
 695                the algorithm in [2]
 696
 697     `Alpha      : Parameter alpha in the formulation (2) or (3). Default is the average of all elements in A. No good justfication for this default value, and you might want to try other values.'
 698
 699     `Beta      : Parameter beta in the formulation (2) or (3).'
 700           Default is the average of all elements in A. No good
 701           justfication for this default value, and you might want to
 702           try other values.
 703
 704     `MaxIter    : Maximum number of iterations. Default is 100.'
 705
 706     `MinIter    : Minimum number of iterations. Default is 20.'
 707
 708     `MaxTime    : Maximum amount of time in seconds. Default is 100,000.'
 709
 710     `Winit      : (m x k) initial value for W.'
 711
 712     `Hinit      : (k x n) initial value for H.'
 713
 714     `Tol        : Stopping tolerance. Default is 1e-3. If you want to obtain a more accurate solution, decrease TOL and increase MAX_ITER at the same time.'
 715
 716     `Verbose  :'
 717
 718          0 (default) - No debugging information is collected.
 719
 720          1 (debugging purpose) - History of computation is returned by 'HIS' variable.
 721
 722          2 (debugging purpose) - History of computation is additionally printed on screen.
 723
 724      *Outputs*
 725     `W    : Obtained basis matrix (m x k)'
 726
 727     `H    : Obtained coefficients matrix (k x n)'
 728
 729     `iter : Number of iterations'
 730
 731     `HIS  : (debugging purpose) History of computation'
 732
 733      Usage Examples:
 734            nmf(A,10)
 735            nmf(A,20,'verbose',2)
 736            nmf(A,30,'verbose',2,'nnls_solver','as')
 737            nmf(A,5,'verbose',2,'type','sparse')
 738            nmf(A,60,'verbose',1,'type','plain','w_init',rand(m,k))
 739            nmf(A,70,'verbose',2,'type','sparse','nnls_solver','bp','alpha',1.1,'beta',1.3)
 740
 741      References:  [1] For using this software, please cite:
 742      Jingu Kim and Haesun Park, Toward Faster Nonnegative Matrix
 743      Factorization: A New Algorithm and Comparisons,
 744      In Proceedings of the 2008 Eighth IEEE International Conference on
 745      Data Mining (ICDM'08), 353-362, 2008
 746      [2] If you use 'nnls_solver'='as' (see below), please cite:
 747      Hyunsoo Kim and Haesun Park, Nonnegative Matrix Factorization
 748      Based
 749      on Alternating Nonnegativity Constrained Least Squares and Active
 750      Set Method,
 751      SIAM Journal on Matrix Analysis and Applications, 2008, 30, 713-730
 752
 753      Check original code at `http://www.cc.gatech.edu/~jingu'
 754
 755      See also: nmf_pg
 756
 757
 758
 759
 760
 761 # name: <cell-element>
 762 # type: sq_string
 763 # elements: 1
 764 # length: 80
 765 Nonnegative Matrix Factorization by Alternating Nonnegativity
 766 Constrained Least
 767
 768
 769
 770 # name: <cell-element>
 771 # type: sq_string
 772 # elements: 1
 773 # length: 6
 774 nmf_pg
 775
 776
 777 # name: <cell-element>
 778 # type: sq_string
 779 # elements: 1
 780 # length: 909
 781  -- Function File: [W, H] = nmf_pg (V, WINIT, HINIT, TOL, TIMELIMIT,
 782           MAXITER)
 783      Non-negative matrix factorization by alternative non-negative
 784      least squares using projected gradients.
 785
 786      The matrix V is factorized into two possitive matrices W and H
 787      such that `V = W*H + U'. Where U is a matrix of residuals that can
 788      be negative or positive. When the matrix V is positive the order
 789      of the elements in U is bounded by the optional named argument TOL
 790      (default value `1e-9').
 791
 792      The factorization is not unique and depends on the inital guess
 793      for the matrices W and H. You can pass this initalizations using
 794      the optional named arguments WINIT and HINIT.
 795
 796      timelimit, maxiter: limit of time and iterations
 797
 798      Examples:
 799
 800             A     = rand(10,5);
 801             [W H] = nmf_pg(A,tol=1e-3);
 802             U     = W*H -A;
 803             disp(max(abs(U)));
 804
 805
 806
 807
 808
 809 # name: <cell-element>
 810 # type: sq_string
 811 # elements: 1
 812 # length: 80
 813 Non-negative matrix factorization by alternative non-negative least
 814 squares usin
 815
 816
 817
 818 # name: <cell-element>
 819 # type: sq_string
 820 # elements: 1
 821 # length: 9
 822 rotparams
 823
 824
 825 # name: <cell-element>
 826 # type: sq_string
 827 # elements: 1
 828 # length: 463
 829  -- Function File: [VSTACKED, ASTACKED] = rotparams (RSTACKED)
 830      The function w = rotparams (r)            - Inverse to rotv().
 831      Using, W    = rotparams(R)  is such that  rotv(w)*r' == eye(3).
 832
 833      If used as, [v,a]=rotparams(r) ,  idem, with v (1 x 3) s.t. w ==
 834      a*v.
 835
 836      0 <= norm(w)==a <= pi
 837
 838      :-O !!  Does not check if 'r' is a rotation matrix.
 839
 840      Ignores matrices with zero rows or with NaNs. (returns 0 for them)
 841
 842      See also: rotv
 843
 844
 845
 846
 847
 848 # name: <cell-element>
 849 # type: sq_string
 850 # elements: 1
 851 # length: 62
 852 The function w = rotparams (r)            - Inverse to rotv().
 853
 854
 855
 856 # name: <cell-element>
 857 # type: sq_string
 858 # elements: 1
 859 # length: 4
 860 rotv
 861
 862
 863 # name: <cell-element>
 864 # type: sq_string
 865 # elements: 1
 866 # length: 500
 867  -- Function File: R =  rotv ( v, ang )
 868      The functionrotv calculates a Matrix of rotation about V w/ angle
 869      |v| r = rotv(v [,ang])
 870
 871      Returns the rotation matrix w/ axis v, and angle, in radians,
 872      norm(v) or ang (if present).
 873
 874      rotv(v) == w'*w + cos(a) * (eye(3)-w'*w) - sin(a) * crossmat(w)
 875
 876      where a = norm (v) and w = v/a.
 877
 878      v and ang may be vertically stacked : If 'v' is 2x3, then rotv( v
 879      ) == [rotv(v(1,:)); rotv(v(2,:))]
 880
 881
 882      See also: rotparams, rota, rot
 883
 884
 885
 886
 887
 888 # name: <cell-element>
 889 # type: sq_string
 890 # elements: 1
 891 # length: 80
 892 The functionrotv calculates a Matrix of rotation about V w/ angle |v| r
 893 = rotv(v
 894
 895
 896
 897 # name: <cell-element>
 898 # type: sq_string
 899 # elements: 1
 900 # length: 8
 901 smwsolve
 902
 903
 904 # name: <cell-element>
 905 # type: sq_string
 906 # elements: 1
 907 # length: 614
 908  -- Function File: X = smwsolve (A, U, V, B)
 909  -- Function File:  smwsolve (SOLVER, U, V, B)
 910      Solves the square system `(A + U*V')*X == B', where U and V are
 911      matrices with several columns, using the Sherman-Morrison-Woodbury
 912      formula, so that a system with A as left-hand side is actually
 913      solved. This is especially advantageous if A is diagonal, sparse,
 914      triangular or positive definite.  A can be sparse or full, the
 915      other matrices are expected to be full.  Instead of a matrix A, a
 916      user may alternatively provide a function SOLVER that performs the
 917      left division operation.
 918
 919
 920
 921
 922 # name: <cell-element>
 923 # type: sq_string
 924 # elements: 1
 925 # length: 80
 926 Solves the square system `(A + U*V')*X == B', where U and V are
 927 matrices with se
 928
 929
 930
 931 # name: <cell-element>
 932 # type: sq_string
 933 # elements: 1
 934 # length: 4
 935 thfm
 936
 937
 938 # name: <cell-element>
 939 # type: sq_string
 940 # elements: 1
 941 # length: 692
 942  -- Function File: Y = thfm (X, MODE)
 943      Trigonometric/hyperbolic functions of square matrix X.
 944
 945      MODE must be the name of a function. Valid functions are 'sin',
 946      'cos', 'tan', 'sec', 'csc', 'cot' and all their inverses and/or
 947      hyperbolic variants, and 'sqrt', 'log' and 'exp'.
 948
 949      The code `thfm (x, 'cos')' calculates matrix cosinus _even if_
 950      input matrix X is _not_ diagonalizable.
 951
 952      _Important note_: This algorithm does _not_ use an eigensystem
 953      similarity transformation. It maps the MODE functions to functions
 954      of `expm', `logm' and `sqrtm', which are known to be robust with
 955      respect to non-diagonalizable ('defective') X.
 956
 957      See also: funm
 958
 959
 960
 961
 962
 963 # name: <cell-element>
 964 # type: sq_string
 965 # elements: 1
 966 # length: 54
 967 Trigonometric/hyperbolic functions of square matrix X.
 968
 969
 970
 971
 972