{"id":6425,"date":"2020-10-23T09:00:51","date_gmt":"2020-10-23T13:00:51","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=6425"},"modified":"2020-12-01T16:14:48","modified_gmt":"2020-12-01T21:14:48","slug":"gil-strang-and-the-cr-matrix-factorization","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2020\/10\/23\/gil-strang-and-the-cr-matrix-factorization\/","title":{"rendered":"Gil Strang and the CR Matrix Factorization"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>My friend Gil Strang is known for his lectures from MIT course 18.06, Linear Algebra, which are available on MIT <a href=\"https:\/\/ocw.mit.edu\/courses\/mathematics\/18-06-linear-algebra-spring-2010\/\">OpenCourseWare<\/a>. He is now describing a new approach to the subject with a series of videos, <a href=\"https:\/\/ocw.mit.edu\/resources\/res-18-010-a-2020-vision-of-linear-algebra-spring-2020\">A 2020 Vision of Linear Algebra<\/a>. This vision is featured in a new book, <a href=\"http:\/\/math.mit.edu\/~gs\/everyone\/\">Linear Algebra for Everyone<\/a>.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#ac2007de-5c98-4a55-8f57-fdc1f061f5ac\">The CR factorization<\/a><\/li><li><a href=\"#4733a7bc-91f3-4c3a-8a51-6b7918e2d566\">Reduced row echelon form<\/a><\/li><li><a href=\"#07c8160b-291a-44ae-9276-9855b68aa706\"><tt>function cr<\/tt><\/a><\/li><li><a href=\"#8745228f-6dc3-45f3-9770-3fbc9c747174\"><tt>magic(4*k)<\/tt><\/a><\/li><li><a href=\"#b7b7ce68-210f-48a8-8c82-481f6cbd3864\"><tt>west0479<\/tt><\/a><\/li><li><a href=\"#bdcebfd5-aba9-4e44-8f41-c938da80083e\">Dueling <tt>rref<\/tt><\/a><\/li><\/ul><\/div><h4>The CR factorization<a name=\"ac2007de-5c98-4a55-8f57-fdc1f061f5ac\"><\/a><\/h4><p>Gil's approach will be familiar to MATLAB users and to readers of this blog.  The key components are matrix factorizations -- LU, QR, eigenvalues and SVD.  But before he gets to those, Gil likes to start with a more fundamental factorization, <tt>A = C*R<\/tt>, that expresses any matrix as a product of a matrix that describes its Column space and a matrix that describes its Row space.<\/p><p>For example, perhaps the first 3-by-3 matrix anybody writes down is<\/p><pre class=\"language-matlab\">A =\r\n<\/pre><pre>      1     2     3\r\n      4     5     6\r\n      7     8     9<\/pre><p>Notice that the middle column is the average of the first and last columns, so the most obvious CR factorization is<\/p><pre class=\"language-matlab\">C =\r\n<\/pre><pre>      1     3\r\n      4     6\r\n      7     9<\/pre><pre class=\"language-matlab\">R =\r\n<\/pre><pre>      1     0.5     0\r\n      0     0.5     1<\/pre><p>However, this not the only possible CR factorization. For another one, keep reading.<\/p><p>The CR factorization works beautifully for the matrices encountered in any introduction to linear algebra.  These matrices are not too large, and their entries are usually small integers.  There are no errors in the input data, and none are expected in the subsequent computation.  But, as Gil freely admits, the CR factorization is not a contender for any serious technical use.<\/p><p>The factorization <tt>A = C*R<\/tt> is <tt>rank_revealing<\/tt>.  The number of columns in <tt>C<\/tt> must be the same as the number of rows in <tt>R<\/tt>. The smallest number of columns for which the product <tt>C*R<\/tt> reproduces <tt>A<\/tt> is <i>defined<\/i> to be the <i>rank<\/i> of <tt>A<\/tt>.  So here, in the first few days of the course, a fundamental concept is introduced.  It goes by many names -- low rank approximation, model reduction, principal components analysis -- and it pervades modern matrix computation. Look under the hood in many AI and ML success stories and you may well find something akin to a CR factorization.<\/p><h4>Reduced row echelon form<a name=\"4733a7bc-91f3-4c3a-8a51-6b7918e2d566\"><\/a><\/h4><p>I think Gil was pleased when he realized that a good choice for <tt>R<\/tt> was an old friend, the <i>reduced row echelon form<\/i>, featured in linear algebra textbooks for decades. Since you've probably forgotten them, or never knew them, these are the requirements for a matrix to be in reduced row echelon form.<\/p><div><ul><li>The first nonzero entry is any nonzero row is a 1, called the leading 1.<\/li><li>All of the other elements in a column containing a leading 1 are zero.<\/li><li>Each leading 1 occurs later in its row than the leading 1's in earlier rows.<\/li><li>Any all zero rows are below the nonzero ones.  (We will soon delete them.)<\/li><\/ul><\/div><p>In the above example, the matrix <tt>R<\/tt> is <i>not<\/i> in reduced row echelon form, because the first nonzero element in the second row is not a 1.<\/p><p>The \"reduced\" in the name refers to the zeros above each leading 1. The algorithm for computing the reduced form is known as Gauss-Jordan elimination.  Its computational complexity is 50 percent more than conventional Gaussian elimination.<\/p><p>An identity matrix is in reduced row echelon form, so if <tt>A<\/tt> is square and invertible, or, in general, has <tt>n<\/tt> linearly independent columns, then <tt>A<\/tt> = <tt>C<\/tt> and <tt>R<\/tt> is the <tt>n<\/tt> -by- <tt>n<\/tt> identity. We expect rank deficient matrices to have more interesting CR factorizations.<\/p><p>One important fact about the reduced form is that it is unique. No matter how you compute it, there is exactly one matrix that meets all the criteria in our bulleted list.  So, if we require <tt>R<\/tt> to be in reduced row echelon form, the CR factorization of any matrix is unique.<\/p><p>You can see the source code for the MATLAB <tt>rref<\/tt> with<\/p><pre class=\"language-matlab\">type <span class=\"string\">rref<\/span>\r\n<\/pre><p>One of the 80 functions in my original pre-MathWorks Fortran MATLAB was <tt>rref<\/tt>, although it didn't come from LINPACK or EISPACK.<\/p><h4><tt>function cr<\/tt><a name=\"07c8160b-291a-44ae-9276-9855b68aa706\"><\/a><\/h4><p>Since MATLAB already has <tt>rref<\/tt>, it took me five minutes, and as many executable lines, to write this function.<\/p><pre class=\"language-matlab\"><span class=\"keyword\">function<\/span> [C,R] = cr(A)\r\n    R = rref(A);\r\n    j = find(sum(double(R)~=0)==1);    <span class=\"comment\">% Indices of leading ones.<\/span>\r\n    r = length(j);                     <span class=\"comment\">% r = rank.<\/span>\r\n    R = R(1:r,:);                      <span class=\"comment\">% Delete all zero rows so R(:,j) == eye(r).<\/span>\r\n    C = A(:,j);\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>Actually, I originally had only four executable lines since <tt>rref<\/tt> also provides the indices <tt>j<\/tt> of the columns containing the leading 1's. But that is only if the class of <tt>A<\/tt> is <tt>single<\/tt> or <tt>double<\/tt>. If <tt>A<\/tt> is a <tt>sym<\/tt>, the Symbolic Math Toolbox makes me compute <tt>j<\/tt> from <tt>R<\/tt>.<\/p><p>The columns of <tt>C<\/tt> are the first <tt>r<\/tt> linearly independent columns of <tt>A<\/tt>. All the linear combination action comes from <tt>R<\/tt>.<\/p><p>The length of <tt>j<\/tt>, the rank, is the most important quantity this function computes.  Everything else depends upon it.<\/p><p>For our example, the above <tt>A<\/tt>,<\/p><pre class=\"language-matlab\">[C,R] = cr(A)\r\n<\/pre><pre class=\"language-matlab\">C =\r\n<\/pre><pre>      1     2\r\n      4     5\r\n      7     8<\/pre><pre class=\"language-matlab\">R =\r\n<\/pre><pre>      1     0    -1\r\n      0     1     2<\/pre><p>Now <tt>C<\/tt> is the first two columns of <tt>A<\/tt> and <tt>R<\/tt> <i>is<\/i> in reduced row echelon form (without the zero row).<\/p><h4><tt>magic(4*k)<\/tt><a name=\"8745228f-6dc3-45f3-9770-3fbc9c747174\"><\/a><\/h4><p>My favorite rank deficient matrices are the MATLAB versions of magic squares with size that is a multiple of four. All of them, no matter how large, have rank 3.  Let's find the CR factorization of <tt>magic(8)<\/tt>.<\/p><pre class=\"codeinput\">   A = magic(8)\r\n   [C,R] = cr(A)\r\n<\/pre><pre class=\"codeoutput\">A =\r\n    64     2     3    61    60     6     7    57\r\n     9    55    54    12    13    51    50    16\r\n    17    47    46    20    21    43    42    24\r\n    40    26    27    37    36    30    31    33\r\n    32    34    35    29    28    38    39    25\r\n    41    23    22    44    45    19    18    48\r\n    49    15    14    52    53    11    10    56\r\n     8    58    59     5     4    62    63     1\r\nC =\r\n    64     2     3\r\n     9    55    54\r\n    17    47    46\r\n    40    26    27\r\n    32    34    35\r\n    41    23    22\r\n    49    15    14\r\n     8    58    59\r\nR =\r\n     1     0     0     1     1     0     0     1\r\n     0     1     0     3     4    -3    -4     7\r\n     0     0     1    -3    -4     4     5    -7\r\n<\/pre><p>The fact that <tt>C<\/tt> has three columns and <tt>R<\/tt> has three rows indicates that the rank of <tt>A<\/tt> is 3. The columns of <tt>C<\/tt> are the first three columns of <tt>A<\/tt>. The first three rows and columns of <tt>R<\/tt> form the 3-by-3 identity matrix. The rest of <tt>R<\/tt> contains the coefficients that generate all of <tt>A<\/tt> from its first three columns.<\/p><p>The patterns that we see in <tt>C<\/tt> and <tt>R<\/tt> come from the algorithm that generates magic squares of order 4*k.  Here are plots of <tt>C<\/tt> and <tt>R<\/tt> for the 24-by-24 magic square.<\/p><pre class=\"language-matlab\">A = magic(24)\r\n[C,R] = cr(A)\r\nplot(C)\r\nplot(R)\r\n<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/magiccr.png\" alt=\"\"> <\/p><h4><tt>west0479<\/tt><a name=\"b7b7ce68-210f-48a8-8c82-481f6cbd3864\"><\/a><\/h4><p>We don't expect the CR factorization to work on larger matrices with inexact entries, but let's try it anyway.<\/p><p>The matrix <tt>west0479<\/tt> has a venerable history. In 1983, Art Westerberg, a chemical engineering professor at Carnegie Mellon, contributed eleven matrices extracted from modelling chemical engineering plants to the Harwell-Boeing Sparse Matrix collection. One of these matrices, <tt>west0479<\/tt>, is a model of an eight-stage distillation column. It is 479-by-479 and has 1887 nonzero entries.<\/p><p>In 1990, when we introduced sparse matrices in MATLAB, we included <tt>west0479<\/tt> in the demos directory.  It has been available in all subsequent releases.<\/p><pre class=\"codeinput\">   clear\r\n   load <span class=\"string\">west0479<\/span>\r\n   whos\r\n<\/pre><pre class=\"codeoutput\">  Name            Size             Bytes  Class     Attributes\r\n\r\n  west0479      479x479            34032  double    sparse    \r\n\r\n<\/pre><p>Today, Google finds 71 web pages that mention <tt>west0479<\/tt>. I stopped counting the number of images on the web that it has generated.<\/p><p>Our <tt>cr<\/tt> will work on the sparse matrix, but it is faster to convert it to full first.<\/p><pre class=\"codeinput\">   A = full(west0479);\r\n   [C,R] = cr(A);\r\n   r = size(C,2)\r\n<\/pre><pre class=\"codeoutput\">r =\r\n   479\r\n<\/pre><p>So, <tt>cr<\/tt> decides this matrix has full rank and returns with <tt>C<\/tt> = <tt>A<\/tt> and <tt>R<\/tt> = <tt>I<\/tt>.  Not very useful.<\/p><h4>Dueling <tt>rref<\/tt><a name=\"bdcebfd5-aba9-4e44-8f41-c938da80083e\"><\/a><\/h4><p>Upon closer examination, we find <tt>rref<\/tt> has an optional second input argument named <tt>tol<\/tt>.  Let's see if that will help us compute a more useful CR factorization.  The <tt>help<\/tt> entry says<\/p><pre>[R,jb] = rref(A,TOL) uses the given tolerance in the rank tests.<\/pre><p>I wrote <tt>rref<\/tt> years ago and I must say now that this description is not very specific.  It should say that any column encountered during the reduction whose largest element is less than or equal to <tt>tol<\/tt> is replaced by all zeros.<\/p><p>The largest element in <tt>west0479<\/tt> is about 3.2e5 and the smallest nonzero is 1.0e-6, so the elements range over eleven orders of magnitude. It turns out that the matrix has many columns with only two nonzero entries, the largest of which is +1 or -1.  So, any <tt>tol<\/tt> larger than 1 causes <tt>rref<\/tt> to ignore those columns.  And any <tt>tol<\/tt> less than 1 includes those columns.  This leads to excessive sparse matrix fill-in.<\/p><p>This animation compares <tt>rref<\/tt> with <tt>tol = 0.99<\/tt> to <tt>rref<\/tt> with <tt>tol = 1.01<\/tt>.  Watch the nonzero counts in the <tt>spy<\/tt> plots.  The count for <tt>tol = 0.99<\/tt> rises to over seven times the count in the original matrix, before it drops somewhat at the end.  The final nonzero count for <tt>tol = 1.01<\/tt> is less than the starting <tt>nnz<\/tt>.  That is not a good sign.<\/p><p>Ultimately, neither value of <tt>tol<\/tt> produces a CR factorization which is close to the original matrix, and no other values do any better. As we expected, CR is just not useful here.<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/west0479_2.gif\" alt=\"\"> <\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_78cb6771038e40fdb44358b2e0059f01() {\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='78cb6771038e40fdb44358b2e0059f01 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 78cb6771038e40fdb44358b2e0059f01';\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 2020 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_78cb6771038e40fdb44358b2e0059f01()\"><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; R2019b<br><\/p><\/div><!--\r\n78cb6771038e40fdb44358b2e0059f01 ##### SOURCE BEGIN #####\r\n%% Gil Strang and the CR Matrix Factorization\r\n% My friend Gil Strang is known for his lectures from MIT course\r\n% 18.06, Linear Algebra, which are available on MIT \r\n% <https:\/\/ocw.mit.edu\/courses\/mathematics\/18-06-linear-algebra-spring-2010\/\r\n% OpenCourseWare>.\r\n% He is now describing a new approach to the subject with a series\r\n% of videos, \r\n% <https:\/\/ocw.mit.edu\/resources\/res-18-010-a-2020-vision-of-linear-algebra-spring-2020\r\n% A 2020 Vision of Linear Algebra>.\r\n% This vision is featured in a new book,\r\n% <http:\/\/math.mit.edu\/~gs\/everyone\/  Linear Algebra for Everyone>.\r\n\r\n%% The CR factorization\r\n% Gil's approach will be familiar to MATLAB users and to readers of\r\n% this blog.  The key components are matrix factorizations REPLACE_WITH_DASH_DASH\r\n% LU, QR, eigenvalues and SVD.  But before he gets to those,\r\n% Gil likes to start with a more fundamental factorization, |A = C*R|, that\r\n% expresses any matrix as a product of a matrix that describes its Column\r\n% space and a matrix that describes its Row space.\r\n%\r\n% For example, perhaps the first 3-by-3 matrix anybody writes down is\r\n%\r\n%   A =\r\n%\r\n%        1     2     3\r\n%        4     5     6\r\n%        7     8     9\r\n% \r\n% Notice that the middle column is the average of the first and last \r\n% columns, so the most obvious CR factorization is\r\n%\r\n%   C =\r\n%\r\n%        1     3\r\n%        4     6\r\n%        7     9\r\n%\r\n%   R =\r\n% \r\n%        1     0.5     0\r\n%        0     0.5     1   \r\n%\r\n% However, this not the only possible CR factorization.\r\n% For another one, keep reading.\r\n% \r\n% The CR factorization works beautifully for the matrices encountered\r\n% in any introduction to linear algebra.  These matrices are not too large,\r\n% and their entries are usually small integers.  There are no errors in\r\n% the input data,\r\n% and none are expected in the subsequent computation.  But, as Gil\r\n% freely admits, the CR factorization is not a contender for any\r\n% serious technical use.\r\n%\r\n% The factorization |A = C*R| is |rank_revealing|.  The\r\n% number of columns in |C| must be the same as the number of rows in |R|.\r\n% The smallest number of columns for which the product |C*R| reproduces\r\n% |A| is _defined_ to be the _rank_ of |A|.  So here, in the first few days\r\n% of the course, a fundamental concept is introduced.  It goes by many \r\n% names REPLACE_WITH_DASH_DASH low rank approximation, model reduction, principal components\r\n% analysis REPLACE_WITH_DASH_DASH and it pervades modern matrix computation. \r\n% Look under the hood in many AI and ML success stories and you may\r\n% well find something akin to a CR factorization.\r\n\r\n%% Reduced row echelon form\r\n% I think Gil was pleased when he realized that a good choice for |R|\r\n% was an old friend, the _reduced row echelon form_, featured in linear\r\n% algebra textbooks for decades.\r\n% Since you've probably forgotten them, or never knew them, these are the\r\n% requirements for a matrix to be in reduced row echelon form.\r\n%\r\n% * The first nonzero entry is any nonzero row is a 1, called the leading 1. \r\n% * All of the other elements in a column containing a leading 1 are zero.\r\n% * Each leading 1 occurs later in its row than the leading 1's in earlier rows.\r\n% * Any all zero rows are below the nonzero ones.  (We will soon delete them.)\r\n%\r\n% In the above example, the matrix |R| is _not_ in reduced row echelon form,\r\n% because the first nonzero element in the second row is not a 1.\r\n%\r\n% The \"reduced\" in the name refers to the zeros above each leading 1.\r\n% The algorithm for computing the reduced form is known as\r\n% Gauss-Jordan elimination.  Its computational complexity is 50 percent\r\n% more than conventional Gaussian elimination.\r\n%\r\n% An identity matrix is in reduced row echelon form, so if |A| is\r\n% square and invertible, or, in general, has |n| linearly independent \r\n% columns, then |A| = |C| and |R| is the |n| -by- |n| identity.\r\n% We expect rank deficient matrices to have more interesting CR factorizations.\r\n%\r\n% One important fact about the reduced form is that it is unique.\r\n% No matter how you compute it, there is exactly one matrix that meets\r\n% all the criteria in our bulleted list.  So, if we require |R| to be\r\n% in reduced row echelon form, the CR factorization of\r\n% any matrix is unique.\r\n%\r\n% You can see the source code for the MATLAB |rref| with\r\n%\r\n%   type rref\r\n% \r\n% One of the 80 functions in my original pre-MathWorks Fortran MATLAB\r\n% was |rref|, although it didn't come from LINPACK or EISPACK.\r\n\r\n%% |function cr|\r\n% Since MATLAB already has |rref|, it took me five minutes, and as many\r\n% executable lines, to write this function.\r\n%\r\n%   function [C,R] = cr(A)\r\n%       R = rref(A);\r\n%       j = find(sum(double(R)~=0)==1);    % Indices of leading ones.\r\n%       r = length(j);                     % r = rank.\r\n%       R = R(1:r,:);                      % Delete all zero rows so R(:,j) == eye(r).\r\n%       C = A(:,j);\r\n%   end\r\n%\r\n% Actually, I originally had only four executable lines since |rref|\r\n% also provides the indices |j| of the columns containing the leading 1's.\r\n% But that is only if the class of |A| is |single| or |double|.\r\n% If |A| is a |sym|, the Symbolic Math Toolbox makes me compute |j|\r\n% from |R|.\r\n%\r\n% The columns of |C| are the first |r| linearly independent columns of |A|.\r\n% All the linear combination action comes from |R|.\r\n%\r\n% The length of |j|, the rank, is the most important quantity\r\n% this function computes.  Everything else depends upon it.\r\n%\r\n% For our example, the above |A|,\r\n%\r\n%   [C,R] = cr(A)\r\n%   \r\n%   C =\r\n%\r\n%        1     2\r\n%        4     5\r\n%        7     8\r\n%\r\n%   R =\r\n%\r\n%        1     0    -1\r\n%        0     1     2\r\n%\r\n% Now |C| is the first two columns of |A| and\r\n% |R| _is_ in reduced row echelon form (without the zero row).\r\n\r\n%% |magic(4*k)|\r\n% My favorite rank deficient matrices are the MATLAB versions\r\n% of magic squares with size that is a multiple of four.\r\n% All of them, no matter how large, have rank 3.  Let's find the CR\r\n% factorization of |magic(8)|.\r\n%\r\n   A = magic(8)\r\n   [C,R] = cr(A) \r\n   \r\n%%\r\n% The fact that |C| has three columns and |R| has three rows\r\n% indicates that the rank of |A| is 3.\r\n% The columns of |C| are the first three columns of |A|.\r\n% The first three rows and columns of |R| form the 3-by-3 identity matrix.\r\n% The rest of |R| contains the coefficients that generate all of |A| from\r\n% its first three columns. \r\n%\r\n% The patterns that we see in |C| and |R| come from the algorithm that\r\n% generates magic squares of order 4*k.  Here are plots of |C| and |R| \r\n% for the 24-by-24 magic square.\r\n%\r\n%   A = magic(24)\r\n%   [C,R] = cr(A)\r\n%   plot(C)\r\n%   plot(R)\r\n%\r\n% <<magiccr.png>>\r\n%\r\n\r\n%% |west0479|\r\n% We don't expect the CR factorization to work on larger matrices with\r\n% inexact entries, but let's try it anyway.\r\n%\r\n% The matrix |west0479| has a venerable history.\r\n% In 1983, Art Westerberg, a chemical engineering professor at\r\n% Carnegie Mellon, contributed eleven matrices extracted from modelling\r\n% chemical engineering plants to the Harwell-Boeing Sparse Matrix\r\n% collection. One of these matrices, |west0479|, is a model of an eight-stage\r\n% distillation column. It is 479-by-479 and has 1887 nonzero entries.\r\n%\r\n% In 1990, when we introduced sparse matrices in MATLAB,\r\n% we included |west0479| in the demos directory.  It has been\r\n% available in all subsequent releases.\r\n\r\n   clear\r\n   load west0479\r\n   whos\r\n\r\n%%\r\n% Today, Google finds 71 web pages that mention |west0479|.\r\n% I stopped counting the number of images on the web that it has generated.\r\n\r\n%%\r\n% Our |cr| will work on the sparse matrix, but it is faster to\r\n% convert it to full first.\r\n\r\n   A = full(west0479);\r\n   [C,R] = cr(A);\r\n   r = size(C,2)\r\n   \r\n%%\r\n% So, |cr| decides this matrix has full rank and returns with |C| = |A|\r\n% and |R| = |I|.  Not very useful.\r\n\r\n%% Dueling |rref|\r\n% Upon closer examination, we find |rref| has an optional second input\r\n% argument named |tol|.  Let's see if that will help us compute\r\n% a more useful CR factorization.  The |help| entry says\r\n%\r\n%  [R,jb] = rref(A,TOL) uses the given tolerance in the rank tests.\r\n% \r\n% I wrote |rref| years ago and I must say now that this description is not\r\n% very specific.  It should say that any column encountered during the\r\n% reduction whose largest element is less than\r\n% or equal to |tol| is replaced by all zeros.\r\n%\r\n% The largest element in |west0479| is about 3.2e5 and the smallest nonzero\r\n% is 1.0e-6, so the elements range over eleven orders of magnitude.\r\n% It turns out that the matrix has many columns with only two\r\n% nonzero entries, the largest of which is +1 or -1.  So, any |tol|\r\n% larger than 1 causes |rref| to ignore those columns.  And any |tol|\r\n% less than 1 includes those columns.  This leads to excessive sparse\r\n% matrix fill-in.\r\n%\r\n% This animation compares |rref| with |tol = 0.99| to |rref| with\r\n% |tol = 1.01|.  Watch the nonzero counts in the |spy| plots.  The count\r\n% for |tol = 0.99| rises to over seven times the count in the original matrix,\r\n% before it drops somewhat at the end.  The final nonzero count for \r\n% |tol = 1.01| is less than the starting |nnz|.  That is not a good sign.\r\n%\r\n% Ultimately,\r\n% neither value of |tol| produces a CR factorization which is close to\r\n% the original matrix, and no other values do any better.\r\n% As we expected, CR is just not useful here.\r\n%\r\n% <<west0479_2.gif>>\r\n%\r\n\r\n##### SOURCE END ##### 78cb6771038e40fdb44358b2e0059f01\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/west0479_small.png\" class=\"img-responsive attachment-post-thumbnail size-post-thumbnail wp-post-image\" alt=\"\" decoding=\"async\" loading=\"lazy\" \/><\/div><!--introduction--><p>My friend Gil Strang is known for his lectures from MIT course 18.06, Linear Algebra, which are available on MIT <a href=\"https:\/\/ocw.mit.edu\/courses\/mathematics\/18-06-linear-algebra-spring-2010\/\">OpenCourseWare<\/a>. He is now describing a new approach to the subject with a series of videos, <a href=\"https:\/\/ocw.mit.edu\/resources\/res-18-010-a-2020-vision-of-linear-algebra-spring-2020\">A 2020 Vision of Linear Algebra<\/a>. This vision is featured in a new book, <a href=\"http:\/\/math.mit.edu\/~gs\/everyone\/\">Linear Algebra for Everyone<\/a>.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2020\/10\/23\/gil-strang-and-the-cr-matrix-factorization\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":6433,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[4,9,6,8],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/6425"}],"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=6425"}],"version-history":[{"count":5,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/6425\/revisions"}],"predecessor-version":[{"id":6512,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/6425\/revisions\/6512"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media\/6433"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=6425"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=6425"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=6425"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}