My friend Gil Strang is known for his lectures from MIT course 18.06, Linear Algebra, which are available on MITOpenCourseWare。He is now describing a new approach to the subject with a series of videos,2020年线性代数的视野。This vision is featured in a new book,Linear Algebra for Everyone。

Contents

这CR factorization

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,A = C*R, that expresses any matrix as a product of a matrix that describes its Column space and a matrix that describes its Row space.

For example, perhaps the first 3-by-3 matrix anybody writes down is

A =

1 2 3 4 5 6 7 8 9

Notice that the middle column is the average of the first and last columns, so the most obvious CR factorization is

C =

1 3 4 6 7 9

r =

1 0.5 0 0 0.5 1

但是，这不是唯一可能的CR分解。对于另一个，请继续阅读。

这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.

这factorizationA = C*Risrank_revealing。这number of columns inCmust be the same as the number of rows inR。产品的最少列数C*RreproducesAisdefined成为rankofA。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.

Reduced row echelon form

I think Gil was pleased when he realized that a good choice forRwas an old friend, thereduced row echelon form, 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.

在上面的示例中，矩阵Risnotin reduced row echelon form, because the first nonzero element in the second row is not a 1.

名称中的“简化”是指每个领先的零1。用于计算还原形式的算法被称为高斯 - 约旦消除。它的计算复杂性比常规高斯消除高50％。

An identity matrix is in reduced row echelon form, so ifAis square and invertible, or, in general, hasn线性独立列，然后A=CandRis then-经过-nidentity. We expect rank deficient matrices to have more interesting CR factorizations.

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 requireR要以降低的行梯形形式，任何矩阵的CR分解都是唯一的。

You can see the source code for the MATLABrrefwith

类型rref

One of the 80 functions in my original pre-MathWorks Fortran MATLAB wasrref，尽管它不是来自Linpack或Eispack。

function cr

由于Matlab已经有rref, it took me five minutes, and as many executable lines, to write this function.

function[C,R] = cr(A) R = rref(A); j = find(sum(double(R)~=0)==1);% Indices of leading ones.r = length(j);％r =等级。r =R(1:r,:);％删除所有零行SO R（：，J）== Eye（r）。C =A(:,j);end

Actually, I originally had only four executable lines sincerref还提供索引j包含领先1的列。但这只有当Aissingleor双倍的。IfAis asym, the Symbolic Math Toolbox makes me computej从R。

这columns ofCare the firstrlinearly independent columns ofA。All the linear combination action comes fromR。

这length ofj, the rank, is the most important quantity this function computes. Everything else depends upon it.

For our example, the aboveA,

[C,R] = cr(A)

C =

1 2 4 5 7 8

r =

1 0 -1 0 1 2

NowCis the first two columns ofAandRis以简化的行梯形形式（无零行）。

magic(4*k)

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魔术（8）。

a =魔术（8）[c，r] = cr（a）

A = 64 2 3 61 60 6 7 57 9 55 54 12 13 51 50 16 17 47 46 20 21 43 42 24 40 26 27 37 36 30 31 33 32 34 35 29 28 38 39 25 41 23 22 44 45 19 18 48 49 15 14 52 53 11 10 56 8 58 59 5 4 62 63 1 C = 64 2 3 9 55 54 17 47 46 40 26 27 32 34 35 41 23 22 49 15 14 8 58 59 R = 1 0 0 1 1 0 0 1 0 1 0 3 4 -3 -4 7 0 0 1 -3 -4 4 5 -7

这fact thatChas three columns andRhas three rows indicates that the rank ofAis 3. The columns ofCare the first three columns ofA。前三行和列Rform the 3-by-3 identity matrix. The rest ofRcontains the coefficients that generate all ofA从its first three columns.

我们看到的模式CandR来自生成4*k的魔法正方形的算法。这是CandRfor the 24-by-24 magic square.

A = magic(24) [C,R] = cr(A) plot(C) plot(R)

west0479

We don't expect the CR factorization to work on larger matrices with inexact entries, but let's try it anyway.

这matrixwest0479有一个古老的历史。1983年,维斯特伯格艺术,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,west0479, is a model of an eight-stage distillation column. It is 479-by-479 and has 1887 nonzero entries.

In 1990, when we introduced sparse matrices in MATLAB, we includedwest0479in the demos directory. It has been available in all subsequent releases.

clear loadwest0479whos

Name Size Bytes Class Attributes west0479 479x479 34032 double sparse

Today, Google finds 71 web pages that mentionwest0479。我停止计算其生成的网络上的图像数量。

我们的crwill work on the sparse matrix, but it is faster to convert it to full first.

A = full(west0479); [C,R] = cr(A); r = size(C,2)

r = 479

So,cr决定此矩阵具有完整的排名，并带有C=AandR=I。Not very useful.

Duelingrref

Upon closer examination, we findrref有一个可选的第二个输入参数tol。让我们看看这是否会帮助我们计算更有用的CR分解。这帮助条目说

[R，JB] = RREF（a，tol）在等级测试中使用给定的公差。

I wroterrefyears 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 totol被所有零取代。

最大的元素west0479is 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, anytol大于1原因rrefto ignore those columns. And anytol少于1包括这些列。这导致过度稀疏的矩阵填充。

This animation comparesrrefwithtol = 0.99torrefwithtol = 1.01。Watch the nonzero counts in thespyplots. The count fortol = 0.99在原始矩阵中的计数上升到末端的七倍以上。最终的非零计数tol = 1.01is less than the startingnnz。That is not a good sign.

最终，两种价值都不tolproduces 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.

Get the MATLAB code

Published with MATLAB® R2019b