{"id":6589,"date":"2021-01-15T07:22:36","date_gmt":"2021-01-15T12:22:36","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=6589"},"modified":"2021-01-15T07:22:36","modified_gmt":"2021-01-15T12:22:36","slug":"cr-and-cab-rank-revealing-matrix-factorizations","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2021\/01\/15\/cr-and-cab-rank-revealing-matrix-factorizations\/","title":{"rendered":"CR and CAB, Rank Revealing Matrix Factorizations"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>The rank of a linear transformation is a fundamental concept in linear algebra and matrix factorizations are fundamental concepts in numerical linear algebra.  Gil Strang's <a href=\"https:\/\/ocw.mit.edu\/resources\/res-18-010-a-2020-vision-of-linear-algebra-spring-2020\/\">2020 Vision of Linear Algebra<\/a> seeks to introduce these notions early in an introductory linear algebra course.<\/p><p>The <tt>CR<\/tt> matrix factorization provides a view of <tt>rref<\/tt>, the reduced row echelon form, as a rank revealing matrix factorization. I discussed <tt>CR<\/tt> in a <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2020\/10\/23\/gil-strang-and-the-cr-matrix-factorization\/\">pair<\/a> of <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2020\/10\/25\/notes-on-cr-and-west0479\/\">posts<\/a> in October. I now want to describe the <tt>CAB<\/tt> factorization, which uses <tt>rref<\/tt> twice in order to treat both rows and columns in the same way.<\/p><p>Gil and I are writing a review paper about these ideas, <i>LU and CR Elimination<\/i>. Here is a link to the current draft, <a href=\"https:\/\/blogs.mathworks.com\/cleve\/files\/LUCR_43.pdf\">LUCR_43.pdf<\/a>.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#000d9b8d-9670-4752-910a-ceb8b2e7a0d9\">CAB<\/a><\/li><li><a href=\"#653966e8-b44f-40de-94fd-b60afe868770\">Rank<\/a><\/li><li><a href=\"#0bdd1f8e-25ea-4e79-9726-dc5599efa29b\">Rank one<\/a><\/li><li><a href=\"#6fb071f8-6f69-4702-873c-8dbbd16ed6e4\">Rank two<\/a><\/li><li><a href=\"#60b6000c-2f55-4930-9e3d-83922c783680\">Full rank<\/a><\/li><li><a href=\"#e75fc345-4b34-4d39-a4d9-8375459efa84\">Franklin magic square<\/a><\/li><li><a href=\"#7b0f9a9f-dd66-47ae-9c90-238a927a8df8\">Backslash and forward slash<\/a><\/li><li><a href=\"#e64af198-f0df-497e-ad5b-7a0d85505db7\">Gauss-Jordan<\/a><\/li><li><a href=\"#558126ca-1170-4f6d-9624-9e70a28d26a1\">Condition<\/a><\/li><li><a href=\"#0606a524-d094-4e36-a15a-4c7d3f90313e\">CAB vs. SVD<\/a><\/li><li><a href=\"#c3b4946b-74c5-4cbb-b38e-251682604e49\">References<\/a><\/li><li><a href=\"#0de5c765-a8b9-4cac-aec0-1d404ef6338e\">Postscript<\/a><\/li><\/ul><\/div><h4>CAB<a name=\"000d9b8d-9670-4752-910a-ceb8b2e7a0d9\"><\/a><\/h4><p>This is a quick description of <tt>cab<\/tt>.<\/p><pre>  [C,W,B] = cab(A) returns C = A(:,cols), B = A(rows,:),\r\n  and W = A(rows,cols) so that<\/pre><pre>      A = C*inv(W)*B<\/pre><pre>  [C,W,B,cols,rows] = cab(A) also returns the indices cols and rows.<\/pre><p>All the elements of <tt>C<\/tt>, <tt>W<\/tt> and <tt>B<\/tt> are taken from the input matrix <tt>A<\/tt> itself; <tt>C<\/tt> is some of its columns, <tt>B<\/tt> is some of its rows, and <tt>W<\/tt> is the set intersection of <tt>C<\/tt> and <tt>B<\/tt>. Actually, <tt>cols<\/tt> and <tt>rows<\/tt> are all this factorization computes. The column indices come from the reduced row echelon form, <tt>rref(A)<\/tt>, and the row indices from the transpose, <tt>rref(A')<\/tt>.<\/p><h4>Rank<a name=\"653966e8-b44f-40de-94fd-b60afe868770\"><\/a><\/h4><p>The fact that the number of linearly independent columns of a matrix is equal to the number of linearly independent rows is not trivial and requires proof.  This number is the <i>rank<\/i> of the matrix and is the size of the nonsingular <tt>W<\/tt>.<\/p><pre>  W is r-by-r where rank(A) = r = length(cols) = length(rows).<\/pre><p>If you try to find <tt>C<\/tt>, <tt>W<\/tt> and <tt>R<\/tt> with a <tt>W<\/tt> that is too small, then <tt>C*inv(W)*R<\/tt> will not be able to reproduce <tt>A<\/tt>.  If <tt>W<\/tt> is too large, it will be singular.<\/p><h4>Rank one<a name=\"0bdd1f8e-25ea-4e79-9726-dc5599efa29b\"><\/a><\/h4><p>If <tt>A<\/tt> has rank one, <tt>A = u*v'<\/tt> for column vector <tt>u<\/tt> and row vector <tt>v'<\/tt>. Then <tt>C<\/tt> is a multiple of <tt>u<\/tt>, <tt>B<\/tt> is a multiple of <tt>v'<\/tt>, and <tt>W<\/tt> is the first nonzero element in <tt>u<\/tt> or <tt>v<\/tt>. Here is an example where the first row is deliberately zero.<\/p><pre class=\"codeinput\">   A = (0:3)'*(4:6)\r\n<\/pre><pre class=\"codeoutput\">A =\r\n     0     0     0\r\n     4     5     6\r\n     8    10    12\r\n    12    15    18\r\n<\/pre><pre class=\"codeinput\">   format <span class=\"string\">compact<\/span>\r\n   [C,W,B] = cab(A)\r\n<\/pre><pre class=\"codeoutput\">C =\r\n     0\r\n     4\r\n     8\r\n    12\r\nW =\r\n     4\r\nB =\r\n     4     5     6\r\n<\/pre><h4>Rank two<a name=\"6fb071f8-6f69-4702-873c-8dbbd16ed6e4\"><\/a><\/h4><p>In this example, even casual inspection will reveal the rank. Note that <tt>A<\/tt> is rectangular.<\/p><pre class=\"codeinput\">   A = reshape(1:24,4,6)'\r\n<\/pre><pre class=\"codeoutput\">A =\r\n     1     2     3     4\r\n     5     6     7     8\r\n     9    10    11    12\r\n    13    14    15    16\r\n    17    18    19    20\r\n    21    22    23    24\r\n<\/pre><pre class=\"codeinput\">   [C,W,B,cols,rows] = cab(A)\r\n<\/pre><pre class=\"codeoutput\">C =\r\n     1     2\r\n     5     6\r\n     9    10\r\n    13    14\r\n    17    18\r\n    21    22\r\nW =\r\n     1     2\r\n     5     6\r\nB =\r\n     1     2     3     4\r\n     5     6     7     8\r\ncols =\r\n     1     2\r\nrows =\r\n     1     2\r\n<\/pre><pre class=\"codeinput\">   rank = length(cols)\r\n<\/pre><pre class=\"codeoutput\">rank =\r\n     2\r\n<\/pre><h4>Full rank<a name=\"60b6000c-2f55-4930-9e3d-83922c783680\"><\/a><\/h4><pre>  If A is square and nonsingular, then C = B = W = A.<\/pre><p>For example,<\/p><pre class=\"codeinput\">   A = pascal(3)\r\n<\/pre><pre class=\"codeoutput\">A =\r\n     1     1     1\r\n     1     2     3\r\n     1     3     6\r\n<\/pre><pre class=\"codeinput\">   [C,W,B] = cab(A)\r\n<\/pre><pre class=\"codeoutput\">C =\r\n     1     1     1\r\n     1     2     3\r\n     1     3     6\r\nW =\r\n     1     1     1\r\n     1     2     3\r\n     1     3     6\r\nB =\r\n     1     1     1\r\n     1     2     3\r\n     1     3     6\r\n<\/pre><h4>Franklin magic square<a name=\"e75fc345-4b34-4d39-a4d9-8375459efa84\"><\/a><\/h4><p>The Franklin <i>semimagic<\/i> square is attributed to Benjamin Franklin in about 1752. A letter from Franklin to one of his colleagues is included in the Franklin papers referenced below.<\/p><p>All the rows and columns of Franklin's square have the required magic sum, but the diagonals do not, so the matrix isn't fully magic. However, many other interesting submatrices are magic, including bent diagonals and any eight elements arranged symmetrically about the center.<\/p><p>What is the rank of this 8-by-8 matrix?<\/p><pre class=\"codeinput\">    A = franklin(8)\r\n<\/pre><pre class=\"codeoutput\">A =\r\n    52    61     4    13    20    29    36    45\r\n    14     3    62    51    46    35    30    19\r\n    53    60     5    12    21    28    37    44\r\n    11     6    59    54    43    38    27    22\r\n    55    58     7    10    23    26    39    42\r\n     9     8    57    56    41    40    25    24\r\n    50    63     2    15    18    31    34    47\r\n    16     1    64    49    48    33    32    17\r\n<\/pre><p>It is hard to think of Franklin's square as a linear transformation, but <tt>cab<\/tt> can still compute its rank.<\/p><pre class=\"codeinput\">   [C,W,B] = cab(A)\r\n<\/pre><pre class=\"codeoutput\">C =\r\n    52    61     4\r\n    14     3    62\r\n    53    60     5\r\n    11     6    59\r\n    55    58     7\r\n     9     8    57\r\n    50    63     2\r\n    16     1    64\r\nW =\r\n    52    61     4\r\n    14     3    62\r\n    53    60     5\r\nB =\r\n    52    61     4    13    20    29    36    45\r\n    14     3    62    51    46    35    30    19\r\n    53    60     5    12    21    28    37    44\r\n<\/pre><p>So, the rank is only three.<\/p><p>Computing <tt>C*inv(W)*B<\/tt> introduces some roundoff error.<\/p><pre class=\"codeinput\">   test = C*inv(W)*B;\r\n   displayer(<span class=\"string\">'test'<\/span>)\r\n<\/pre><pre class=\"codeoutput\">test =\r\n    52.00  61.00   4.00  13.00  20.00  29.00  36.00  45.00\r\n    14.00   3.00  62.00  51.00  46.00  35.00  30.00  19.00\r\n    53.00  60.00   5.00  12.00  21.00  28.00  37.00  44.00\r\n    11.00   6.00  59.00  54.00  43.00  38.00  27.00  22.00\r\n    55.00  58.00   7.00  10.00  23.00  26.00  39.00  42.00\r\n     9.00   8.00  57.00  56.00  41.00  40.00  25.00  24.00\r\n    50.00  63.00   2.00  15.00  18.00  31.00  34.00  47.00\r\n    16.00   1.00  64.00  49.00  48.00  33.00  32.00  17.00\r\n<\/pre><p>Rounding <tt>test<\/tt> to the nearest integers reproduces Franklin's magic square exactly.<\/p><h4>Backslash and forward slash<a name=\"7b0f9a9f-dd66-47ae-9c90-238a927a8df8\"><\/a><\/h4><pre>  Computing C*inv(W)*B by (C\/W)*B or C*(W\\B) is usually more accurate.\r\n  If A has integer elements, you might try round(C\/W*B) or round(C*(W\\B)).<\/pre><p>This 5-by-5 matrix was deliberately constructed so that its rank deficiency is not obvious.<\/p><pre class=\"codeinput\">   A  = gallery(5)\r\n<\/pre><pre class=\"codeoutput\">A =\r\n          -9          11         -21          63        -252\r\n          70         -69         141        -421        1684\r\n        -575         575       -1149        3451      -13801\r\n        3891       -3891        7782      -23345       93365\r\n        1024       -1024        2048       -6144       24572\r\n<\/pre><pre class=\"codeinput\">   [C,W,B] = cab(A)\r\n<\/pre><pre class=\"codeoutput\">C =\r\n          -9          11         -21          63\r\n          70         -69         141        -421\r\n        -575         575       -1149        3451\r\n        3891       -3891        7782      -23345\r\n        1024       -1024        2048       -6144\r\nW =\r\n          -9          11         -21          63\r\n          70         -69         141        -421\r\n        -575         575       -1149        3451\r\n        3891       -3891        7782      -23345\r\nB =\r\n          -9          11         -21          63        -252\r\n          70         -69         141        -421        1684\r\n        -575         575       -1149        3451      -13801\r\n        3891       -3891        7782      -23345       93365\r\n<\/pre><pre class=\"codeinput\">   rank = size(W,1)\r\n<\/pre><pre class=\"codeoutput\">rank =\r\n     4\r\n<\/pre><p>In this example, <tt>W<\/tt> is mildly ill-conditioned.<\/p><pre class=\"codeinput\">   kappa = cond(W)\r\n<\/pre><pre class=\"codeoutput\">kappa =\r\n   2.4301e+04\r\n<\/pre><p>There are three ways to recreate <tt>A<\/tt>.<\/p><pre class=\"codeinput\">   invw = C*inv(W)*B;\r\n   back = (C\/W) * B;\r\n   fore = C * (W\\B);\r\n<\/pre><p>The relative error with inversion is affected by the condition of <tt>W<\/tt>, while backslash and forward slash are almost five orders of magnitude more accurate,<\/p><pre class=\"codeinput\">   format <span class=\"string\">short<\/span> <span class=\"string\">e<\/span>\r\n   relerr = [norm(A-invw) norm(A-fore) norm(A-back)]\/norm(A)\r\n<\/pre><pre class=\"codeoutput\">relerr =\r\n   2.1266e-13   4.7760e-18   3.7217e-17\r\n<\/pre><p>Rounding any of the three to integers reproduces <tt>A<\/tt> exactly.<\/p><pre class=\"codeinput\">   A = round(C*inv(W)*B)\r\n<\/pre><pre class=\"codeoutput\">A =\r\n          -9          11         -21          63        -252\r\n          70         -69         141        -421        1684\r\n        -575         575       -1149        3451      -13801\r\n        3891       -3891        7782      -23345       93365\r\n        1024       -1024        2048       -6144       24572\r\n<\/pre><h4>Gauss-Jordan<a name=\"e64af198-f0df-497e-ad5b-7a0d85505db7\"><\/a><\/h4><p>The algorithm for computing the <tt>rref<\/tt> is Gauss-Jordan elimination, which has a dubious reputation in the numerical analysis community.  After listening to me trying to characterize the algorithm's notoriety, Gil wrote in the introductory paragraphs of <a href=\"https:\/\/blogs.mathworks.com\/cleve\/files\/LUCR_43.pdf\">our survey<\/a> that Gauss-Jordan is \"expensive and numerically perilous\".<\/p><p>The \"expensive\" part is easy to explain.  The computation of the <tt>rref<\/tt> introduces zeros in a column both above and below the pivot element, while Gaussian elimination only zeros below the pivot. That's less work. We can measure computational complexity in terms of a fused multiply-add instruction, FMA, that multiplies two floating point values and then adds a third value to the result in one operation. For a large n-by-n system, Gauss-Jordan requires about (1\/2) n^3 FMAs, while Gaussian elimination requires less, only (1\/3) n^3 FMAs.<\/p><p>However, we don't intend <tt>CR<\/tt> and <tt>CAB<\/tt> to be used with large matrices. For our purposes, even using Gauss-Jordan twice is not expensive.<\/p><p>The \"numerically perilous\" portion of Gauss-Jordan's reputation is more subtle.  It derives from solving a linear system <tt>A*x = b<\/tt> by computing the <tt>rref<\/tt> of an augmented matrix, <tt>[A b]<\/tt>. Gaussian elimination, on the other hand, finds triangular factors <tt>L<\/tt> and <tt>U<\/tt> and then computes <tt>x<\/tt> by back substitution. We could say that Gaussian elimination is <i>partial<\/i> elimination, while Gauss-Jordan is <i>complete<\/i> elimination.  The two approaches are compared in the 1975 paper by Jim Wilkinson and Gwen Peters referenced below. Gaussian elimination almost always computes the solution of a nearby system, but the same cannot be said for Gauss-Jordan on the augmented matrix.<\/p><p>However, we are not using either form of elimination to actually compute the <tt>CAB<\/tt> factors; <tt>rref<\/tt> just selects indices.  The matrices <tt>C<\/tt>, <tt>W<\/tt> and <tt>B<\/tt> from <tt>CAB(A)<\/tt> are submatrices of <tt>A<\/tt>.  They are certainly accurate.  There are two potential numerical difficulties -- <tt>W<\/tt> might be badly conditioned or it might be the wrong size.<\/p><h4>Condition<a name=\"558126ca-1170-4f6d-9624-9e70a28d26a1\"><\/a><\/h4><p>What if <tt>A<\/tt> is badly conditioned?  Then <tt>W<\/tt> has to inherit that conditioning.  Here is an example.  Start with<\/p><pre class=\"codeinput\">   n = 12;\r\n   U = eye(n,n) - triu(ones(n,n),1)\r\n<\/pre><pre class=\"codeoutput\">U =\r\n     1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1\r\n     0     1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1\r\n     0     0     1    -1    -1    -1    -1    -1    -1    -1    -1    -1\r\n     0     0     0     1    -1    -1    -1    -1    -1    -1    -1    -1\r\n     0     0     0     0     1    -1    -1    -1    -1    -1    -1    -1\r\n     0     0     0     0     0     1    -1    -1    -1    -1    -1    -1\r\n     0     0     0     0     0     0     1    -1    -1    -1    -1    -1\r\n     0     0     0     0     0     0     0     1    -1    -1    -1    -1\r\n     0     0     0     0     0     0     0     0     1    -1    -1    -1\r\n     0     0     0     0     0     0     0     0     0     1    -1    -1\r\n     0     0     0     0     0     0     0     0     0     0     1    -1\r\n     0     0     0     0     0     0     0     0     0     0     0     1\r\n<\/pre><p><tt>U<\/tt> is badly conditioned. <tt>cond(U,1)<\/tt> and <tt>cond(U,inf)<\/tt> are both equal to <tt>n*2^(n-1)<\/tt>.<\/p><p>Let<\/p><pre class=\"codeinput\">   A = U'*U\r\n<\/pre><pre class=\"codeoutput\">A =\r\n     1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1\r\n    -1     2     0     0     0     0     0     0     0     0     0     0\r\n    -1     0     3     1     1     1     1     1     1     1     1     1\r\n    -1     0     1     4     2     2     2     2     2     2     2     2\r\n    -1     0     1     2     5     3     3     3     3     3     3     3\r\n    -1     0     1     2     3     6     4     4     4     4     4     4\r\n    -1     0     1     2     3     4     7     5     5     5     5     5\r\n    -1     0     1     2     3     4     5     8     6     6     6     6\r\n    -1     0     1     2     3     4     5     6     9     7     7     7\r\n    -1     0     1     2     3     4     5     6     7    10     8     8\r\n    -1     0     1     2     3     4     5     6     7     8    11     9\r\n    -1     0     1     2     3     4     5     6     7     8     9    12\r\n<\/pre><p><tt>A<\/tt> is symmetric, positive definite and full rank. The determinant of <tt>A<\/tt> is equal to one.  I don't see a neat formula for <tt>cond(A)<\/tt>, but I know it does grow like <tt>4^n<\/tt>.<\/p><p>Since <tt>A<\/tt> is full rank, its <tt>A = C*inv(W)*B<\/tt> factorization must be <tt>A = A*inv(A)*A<\/tt>.<\/p><pre class=\"codeinput\">   [C,W,B] = cab(A);\r\n   test = [isequal(C,A) isequal(W,A) isequal(B,A)]\r\n<\/pre><pre class=\"codeoutput\">test =\r\n  1&times;3 logical array\r\n   1   1   1\r\n<\/pre><p><tt>W<\/tt> is badly conditioned.  <tt>inv(W)<\/tt> has large elements.<\/p><pre class=\"codeinput\">   Winv = inv(W)\r\n<\/pre><pre class=\"codeoutput\">Winv =\r\n  Columns 1 through 6\r\n     1398102      699051      349526      174764       87384       43696\r\n      699051      349526      174763       87382       43692       21848\r\n      349526      174763       87382       43691       21846       10924\r\n      174764       87382       43691       21846       10923        5462\r\n       87384       43692       21846       10923        5462        2731\r\n       43696       21848       10924        5462        2731        1366\r\n       21856       10928        5464        2732        1366         683\r\n       10944        5472        2736        1368         684         342\r\n        5504        2752        1376         688         344         172\r\n        2816        1408         704         352         176          88\r\n        1536         768         384         192          96          48\r\n        1024         512         256         128          64          32\r\n  Columns 7 through 12\r\n       21856       10944        5504        2816        1536        1024\r\n       10928        5472        2752        1408         768         512\r\n        5464        2736        1376         704         384         256\r\n        2732        1368         688         352         192         128\r\n        1366         684         344         176          96          64\r\n         683         342         172          88          48          32\r\n         342         171          86          44          24          16\r\n         171          86          43          22          12           8\r\n          86          43          22          11           6           4\r\n          44          22          11           6           3           2\r\n          24          12           6           3           2           1\r\n          16           8           4           2           1           1\r\n<\/pre><pre class=\"codeinput\">   logWinv = floor(log2(Winv))\r\n<\/pre><pre class=\"codeoutput\">logWinv =\r\n    20    19    18    17    16    15    14    13    12    11    10    10\r\n    19    18    17    16    15    14    13    12    11    10     9     9\r\n    18    17    16    15    14    13    12    11    10     9     8     8\r\n    17    16    15    14    13    12    11    10     9     8     7     7\r\n    16    15    14    13    12    11    10     9     8     7     6     6\r\n    15    14    13    12    11    10     9     8     7     6     5     5\r\n    14    13    12    11    10     9     8     7     6     5     4     4\r\n    13    12    11    10     9     8     7     6     5     4     3     3\r\n    12    11    10     9     8     7     6     5     4     3     2     2\r\n    11    10     9     8     7     6     5     4     3     2     1     1\r\n    10     9     8     7     6     5     4     3     2     1     1     0\r\n    10     9     8     7     6     5     4     3     2     1     0     0\r\n<\/pre><p>Follow the diagonal of <tt>logWinv<\/tt> from <tt>(n,n)<\/tt> to <tt>(1,1)<\/tt>.<\/p><pre class=\"codeinput\">   Winv(1,1) &gt; 4^(n-2)\r\n<\/pre><pre class=\"codeoutput\">ans =\r\n  logical\r\n   1\r\n<\/pre><p>Is <tt>A<\/tt> close to a matrix of lower rank?  The answer is yes, but the lower rank approximation comes from the SVD.<\/p><h4>CAB vs. SVD<a name=\"0606a524-d094-4e36-a15a-4c7d3f90313e\"><\/a><\/h4><p>The <tt>SVD<\/tt> is the ultimate rank revealing factorization, but it does not actually attempt to compute the rank.<\/p><pre>  [U,S,V] = svd(A)<\/pre><p>computes a full set of singular values and vectors as if the matrix had full rank. The singular values, <tt>s = diag(S)<\/tt>, are indexed in decreasing order,<\/p><pre class=\"language-matlab\">s(1) &gt;=  . . . &gt;= s(r) &gt;= s(r+1) . . .\r\n<\/pre><p>You can compare the singular values with a tolerance to obtain an approximate rank and factorization.<\/p><p>The choice of an effective rank is backed by the <a href=\"https:\/\/en.wikipedia.org\/wiki\/Low-rank_approximation\">Eckart-Young-Mirsky<\/a> theorem: If <tt>Ar<\/tt> is the approximation made by choosing the rank to be <tt>r<\/tt>,<\/p><pre class=\"language-matlab\">Ar = U(:,1:r)*S(1:r,1:r)*V(:,1:r)'\r\n<\/pre><p>the norm of the error is bounded by the first term that is omitted,<\/p><pre class=\"language-matlab\">norm(A - Ar,2) &lt;= s(r+1)\r\n<\/pre><p>If the singular values decrease rapidly, an accurate low rank approximation is achievable with only a few terms.<\/p><p>On the other hand, <tt>CAB<\/tt> chooses bases from the columns and rows of a matrix itself. There is no Eckert-Young theorem for the <tt>CR<\/tt> or <tt>CAB<\/tt> factorizations.<\/p><p>We intend <tt>CR<\/tt> and <tt>CAB<\/tt> to be used in the classroom on small matrices, usually with integer entries.  We will let the SVD and randomized matrix factorizations handle the big stuff.<\/p><h4>References<a name=\"c3b4946b-74c5-4cbb-b38e-251682604e49\"><\/a><\/h4><p>Benjamin Franklin, <i>The Papers of Benjamin Franklin<\/i>, Volume 44, 1750-53, To Peter Collinson. (American Philosophical Society and Yale University) p. 392. <a href=\"https:\/\/franklinpapers.org\/framedVolumes.jsp?vol=4&amp;page=392a\">https:\/\/franklinpapers.org\/framedVolumes.jsp?vol=4&amp;page=392a<\/a>.<\/p><p>G. Peters and J. H. Wilkinson, \"On the stability of Gauss-Jordan elimination with pivoting\", <i>Communications ACM<\/i> 18, 1975, p. 20-24. <a href=\"https:\/\/dl.acm.org\/doi\/10.1145\/360569.360653\">https:\/\/dl.acm.org\/doi\/10.1145\/360569.360653<\/a>, also at <a href=\"http:\/\/www.stat.uchicago.edu\/~lekheng\/courses\/302\/gauss-jordan2.pdf\">http:\/\/www.stat.uchicago.edu\/~lekheng\/courses\/302\/gauss-jordan2.pdf<\/a><\/p><h4>Postscript<a name=\"0de5c765-a8b9-4cac-aec0-1d404ef6338e\"><\/a><\/h4><p>Ben Franklin can also supply the graphic for today's post.<\/p><pre class=\"codeinput\">    A = franklin(8);\r\n    [~,~,B] = cab(A);\r\n    plot(B)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/cab_blog_01.png\" alt=\"\"> <script language=\"JavaScript\"> <!-- \r\n    function grabCode_409a020ca4d84f25915a1dbd0cf3358d() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='409a020ca4d84f25915a1dbd0cf3358d ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 409a020ca4d84f25915a1dbd0cf3358d';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2021 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_409a020ca4d84f25915a1dbd0cf3358d()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2021a<br><\/p><\/div><!--\r\n409a020ca4d84f25915a1dbd0cf3358d ##### SOURCE BEGIN #####\r\n%% CR and CAB, Rank Revealing Matrix Factorizations\r\n% The rank of a linear transformation is a fundamental concept in linear\r\n% algebra and matrix factorizations are fundamental concepts in numerical\r\n% linear algebra.  Gil Strang's\r\n% <https:\/\/ocw.mit.edu\/resources\/res-18-010-a-2020-vision-of-linear-algebra-spring-2020\/\r\n% 2020 Vision of Linear Algebra> seeks to introduce these notions early\r\n% in an introductory linear algebra course.\r\n%\r\n% The |CR| matrix factorization provides a view of |rref|, the \r\n% reduced row echelon form, as a rank revealing matrix factorization.\r\n% I discussed |CR| in a\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2020\/10\/23\/gil-strang-and-the-cr-matrix-factorization\/\r\n% pair> of\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2020\/10\/25\/notes-on-cr-and-west0479\/\r\n% posts> in October.\r\n% I now want to describe the |CAB| factorization, which uses |rref| twice\r\n% in order to treat both rows and columns in the same way.\r\n%\r\n% Gil and I are writing a review paper about these ideas,\r\n% _LU and CR Elimination_.\r\n% Here is a link to the current draft, \r\n% <https:\/\/blogs.mathworks.com\/cleve\/files\/LUCR_43.pdf\r\n% LUCR_43.pdf>.\r\n\r\n%% CAB\r\n% This is a quick description of |cab|.\r\n%\r\n%    [C,W,B] = cab(A) returns C = A(:,cols), B = A(rows,:),\r\n%    and W = A(rows,cols) so that\r\n%\r\n%        A = C*inv(W)*B\r\n%\r\n%    [C,W,B,cols,rows] = cab(A) also returns the indices cols and rows.\r\n%\r\n% All the elements of |C|, |W| and |B| are taken from the input matrix\r\n% |A| itself; |C| is some of its columns, |B| is some of its rows, \r\n% and |W| is the set intersection of |C| and |B|.\r\n% Actually, |cols| and |rows| are all this factorization computes.\r\n% The column indices come from the reduced row echelon form, |rref(A)|,\r\n% and the row indices from the transpose, |rref(A')|.\r\n\r\n%% Rank\r\n% The fact that the number of linearly independent columns of a matrix\r\n% is equal to the number of linearly independent rows is not trivial and\r\n% requires proof.  This number is the _rank_ of the matrix and is the\r\n% size of the nonsingular |W|.\r\n%   \r\n%    W is r-by-r where rank(A) = r = length(cols) = length(rows).\r\n%\r\n% If you try to find |C|, |W| and |R| with a |W| that is too small,\r\n% then |C*inv(W)*R| will not be able to reproduce |A|.  If |W|\r\n% is too large, it will be singular.\r\n\r\n%% Rank one\r\n% If |A| has rank one, |A = u*v'| for column vector |u| and row vector |v'|.\r\n% Then |C| is a multiple of |u|, |B| is a multiple of |v'|,\r\n% and |W| is the first nonzero element in |u| or |v|.\r\n% Here is an example where the first row is deliberately zero.\r\n%\r\n   A = (0:3)'*(4:6)\r\n   \r\n%%\r\n   format compact\r\n   [C,W,B] = cab(A)\r\n\r\n%% Rank two\r\n% In this example, even casual inspection will reveal the rank.\r\n% Note that |A| is rectangular.\r\n\r\n   A = reshape(1:24,4,6)'\r\n   \r\n%%\r\n   \r\n   [C,W,B,cols,rows] = cab(A)\r\n\r\n%%\r\n   \r\n   rank = length(cols)\r\n   \r\n%% Full rank\r\n%\r\n%    If A is square and nonsingular, then C = B = W = A.\r\n%\r\n% For example,\r\n% \r\n   A = pascal(3)\r\n   \r\n%%\r\n   [C,W,B] = cab(A)   \r\n\r\n%% Franklin magic square\r\n% The Franklin _semimagic_ square is attributed to Benjamin Franklin\r\n% in about 1752. A letter from Franklin to one of his colleagues\r\n% is included in the Franklin papers referenced below.\r\n%\r\n% All the rows and columns of Franklin's square have the required magic\r\n% sum, but the diagonals do not, so the matrix isn't fully magic.\r\n% However, many other interesting submatrices are magic, including bent\r\n% diagonals and any eight elements arranged symmetrically about the center.\r\n%\r\n% What is the rank of this 8-by-8 matrix?\r\n\r\n    A = franklin(8)\r\n    \r\n%%\r\n% It is hard to think of Franklin's square as a linear transformation,\r\n% but |cab| can still compute its rank.\r\n\r\n   [C,W,B] = cab(A)\r\n\r\n%%\r\n% So, the rank is only three.\r\n%\r\n% Computing |C*inv(W)*B| introduces some roundoff error.\r\n\r\n   test = C*inv(W)*B;\r\n   displayer('test')\r\n   \r\n%%\r\n% Rounding |test| to the nearest integers reproduces Franklin's magic\r\n% square exactly.\r\n\r\n%% Backslash and forward slash\r\n%\r\n%    Computing C*inv(W)*B by (C\/W)*B or C*(W\\B) is usually more accurate.\r\n%    If A has integer elements, you might try round(C\/W*B) or round(C*(W\\B)).\r\n%\r\n% This 5-by-5 matrix was deliberately constructed so that its rank deficiency\r\n% is not obvious.\r\n\r\n   A  = gallery(5)\r\n   \r\n%%\r\n   [C,W,B] = cab(A)\r\n   \r\n%%\r\n   rank = size(W,1)\r\n\r\n%%\r\n% In this example, |W| is mildly ill-conditioned.\r\n\r\n   kappa = cond(W)\r\n   \r\n%%\r\n% There are three ways to recreate |A|.\r\n\r\n   invw = C*inv(W)*B;\r\n   back = (C\/W) * B;\r\n   fore = C * (W\\B);\r\n   \r\n%%\r\n% The relative error with inversion is affected by the condition of |W|,\r\n% while backslash and forward slash are almost five orders of magnitude\r\n% more accurate,\r\n\r\n   format short e\r\n   relerr = [norm(A-invw) norm(A-fore) norm(A-back)]\/norm(A)\r\n   \r\n%%\r\n% Rounding any of the three to integers reproduces |A| exactly.\r\n\r\n   A = round(C*inv(W)*B)\r\n   \r\n%% Gauss-Jordan\r\n% The algorithm for computing the |rref| is Gauss-Jordan elimination,\r\n% which has a dubious reputation in the numerical analysis\r\n% community.  After listening to me trying to characterize the\r\n% algorithm's notoriety, Gil wrote in the introductory paragraphs of\r\n% <https:\/\/blogs.mathworks.com\/cleve\/files\/LUCR_43.pdf\r\n% our survey> that Gauss-Jordan is \"expensive and numerically perilous\".\r\n%\r\n% The \"expensive\" part is easy to explain.  The computation of the |rref| \r\n% introduces zeros in a column both above and below the pivot element,\r\n% while Gaussian elimination only zeros below the pivot.\r\n% That's less work.\r\n% We can measure computational complexity in terms of a\r\n% fused multiply-add instruction, FMA,\r\n% that multiplies two floating point values and then adds a third value \r\n% to the result in one operation. For a large n-by-n system, Gauss-Jordan \r\n% requires about (1\/2) n^3 FMAs, while Gaussian elimination requires less,\r\n% only (1\/3) n^3 FMAs. \r\n%\r\n% However, we don't intend |CR| and |CAB| to be used with large matrices.\r\n% For our purposes, even using Gauss-Jordan twice is not expensive.\r\n%\r\n% The \"numerically perilous\" portion of Gauss-Jordan's reputation\r\n% is more subtle.  It derives from solving a linear system |A*x = b| \r\n% by computing the |rref| of an augmented matrix, |[A b]|.\r\n% Gaussian elimination, on the other hand, finds triangular factors\r\n% |L| and |U| and then computes |x| by back substitution.\r\n% We could say that Gaussian elimination is _partial_ elimination,\r\n% while Gauss-Jordan is _complete_ elimination.  The two\r\n% approaches are compared in the 1975 paper by Jim Wilkinson and Gwen Peters\r\n% referenced below. Gaussian elimination almost always computes the solution\r\n% of a nearby system, but the same cannot be said for Gauss-Jordan on the\r\n% augmented matrix.\r\n%\r\n% However, we are not using either form of elimination to actually compute\r\n% the |CAB| factors; |rref| just selects indices.  The matrices |C|, |W|\r\n% and |B| from |CAB(A)| are submatrices of |A|.  They are certainly \r\n% accurate.  There are two potential numerical difficulties REPLACE_WITH_DASH_DASH |W| might\r\n% be badly conditioned or it might be the wrong size.\r\n\r\n%% Condition\r\n% What if |A| is badly conditioned?  Then |W| has to inherit that\r\n% conditioning.  Here is an example.  Start with\r\n\r\n   n = 12;\r\n   U = eye(n,n) - triu(ones(n,n),1)\r\n   \r\n%%\r\n% |U| is badly conditioned. |cond(U,1)| and |cond(U,inf)| are\r\n% both equal to |n*2^(n-1)|.\r\n   \r\n%%\r\n% Let\r\n\r\n   A = U'*U\r\n   \r\n%%\r\n% |A| is symmetric, positive definite and full rank.\r\n% The determinant of |A| is equal to one.  I don't see a neat formula\r\n% for |cond(A)|, but I know it does grow like |4^n|.\r\n   \r\n%%\r\n% Since |A| is full rank, its |A = C*inv(W)*B| factorization must be\r\n% |A = A*inv(A)*A|.  \r\n\r\n   [C,W,B] = cab(A);\r\n   test = [isequal(C,A) isequal(W,A) isequal(B,A)]\r\n   \r\n%%\r\n% |W| is badly conditioned.  |inv(W)| has large elements.\r\n\r\n   Winv = inv(W)\r\n   \r\n%%\r\n\r\n   logWinv = floor(log2(Winv))\r\n   \r\n%%\r\n% Follow the diagonal of |logWinv| from |(n,n)| to |(1,1)|.\r\n\r\n   Winv(1,1) > 4^(n-2)\r\n   \r\n%%\r\n% Is |A| close to a matrix of lower rank?  The answer is yes,\r\n% but the lower rank approximation comes from the SVD.\r\n\r\n%% CAB vs. SVD\r\n% The |SVD| is the ultimate rank revealing factorization, but it does not\r\n% actually attempt to compute the rank.\r\n%\r\n%    [U,S,V] = svd(A)\r\n%\r\n% computes a full set of singular values and vectors\r\n% as if the matrix had full rank.\r\n% The singular values, |s = diag(S)|, are indexed in decreasing order,\r\n%\r\n%   s(1) >=  . . . >= s(r) >= s(r+1) . . . \r\n%\r\n% You can compare the singular values with a tolerance to\r\n% obtain an approximate rank and factorization.\r\n\r\n\r\n%%\r\n% The choice of an effective rank is backed by the\r\n% <https:\/\/en.wikipedia.org\/wiki\/Low-rank_approximation\r\n% Eckart-Young-Mirsky> theorem:\r\n% If |Ar| is the approximation made by choosing the rank to be |r|,\r\n%\r\n%   Ar = U(:,1:r)*S(1:r,1:r)*V(:,1:r)'\r\n%\r\n% the norm of the error is bounded by the first term that is\r\n% omitted,\r\n%\r\n%   norm(A - Ar,2) <= s(r+1)\r\n%\r\n% If the singular values decrease rapidly, an accurate low rank\r\n% approximation is achievable with only a few terms.\r\n\r\n\r\n%%\r\n% On the other hand, |CAB| chooses bases from the columns and rows\r\n% of a matrix itself.  \r\n% There is no Eckert-Young theorem for the |CR| or |CAB| factorizations.\r\n%\r\n% We intend |CR| and |CAB| to be used in the classroom on small matrices,\r\n% usually with integer entries.  We will let the SVD and randomized matrix\r\n% factorizations handle the big stuff.\r\n\r\n%% References\r\n% Benjamin Franklin, _The Papers of Benjamin Franklin_,\r\n% Volume 44, 1750-53, To Peter Collinson.\r\n% (American Philosophical Society and Yale University) p. 392.\r\n% <https:\/\/franklinpapers.org\/framedVolumes.jsp?vol=4&page=392a>.\r\n%\r\n% G. Peters and J. H. Wilkinson,\r\n% \"On the stability of Gauss-Jordan elimination with pivoting\",\r\n% _Communications ACM_ 18, 1975, p. 20-24.\r\n% <https:\/\/dl.acm.org\/doi\/10.1145\/360569.360653>,\r\n% also at\r\n% <http:\/\/www.stat.uchicago.edu\/~lekheng\/courses\/302\/gauss-jordan2.pdf>\r\n\r\n%% Postscript\r\n% Ben Franklin can also supply the graphic for today's post.\r\n\r\n    A = franklin(8);\r\n    [~,~,B] = cab(A);\r\n    plot(B)\r\n    \r\n\r\n##### SOURCE END ##### 409a020ca4d84f25915a1dbd0cf3358d\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/cab_blog_01.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>The rank of a linear transformation is a fundamental concept in linear algebra and matrix factorizations are fundamental concepts in numerical linear algebra.  Gil Strang's <a href=\"https:\/\/ocw.mit.edu\/resources\/res-18-010-a-2020-vision-of-linear-algebra-spring-2020\/\">2020 Vision of Linear Algebra<\/a> seeks to introduce these notions early in an introductory linear algebra course.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2021\/01\/15\/cr-and-cab-rank-revealing-matrix-factorizations\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[9,6,16,30],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/6589"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=6589"}],"version-history":[{"count":1,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/6589\/revisions"}],"predecessor-version":[{"id":6591,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/6589\/revisions\/6591"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=6589"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=6589"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=6589"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}