What is A\A?

The answer:A\Ais alwaysI, except when it isn't.

Contents

Why A\A ?

I have been explaining our backslash operator for almost 50 years, but I have to admit that the title of today's blog post looks a bit strange. You never see3\3for numbers. So, what isA\A?

A\Asolves the equation

|A*X = A|

IfAis square and nonsingular, thenX=I。But what ifAis rectangular, or singular, or not known exactly? These are all nontrivial questions.

This post is my response to a recent internal discussion at MathWorks about backslash generating NaNs.

Mono-elemental matrices

Any general statement about matrices should be applicable to 1-by-1 matrices in particular. For 1-by-1 matrices, there is an easy answer to our question.

  • Ifais any nonzero number, thena\a = 1

My colleague Pete Stewart likes to use "mono-elemental" for this important class of matrices.

Nonsingular matrices

When the 1-by-1 case is generalized to n-by-n with larger n, it becomes:

  • IfAis any nonsingular matrix, thenA\A = I,

This is not the end of the story, of course. It's our job to investigate approximately nonsingular and approximately equals.

Mono-elemental再次

When my daughter was in the fifth grade, her math teacher told her that mathematicians hadn't figured out yet how to divide by zero. But the authors of the IEEE 754 standard for floating point arithmetic have figured it out and have assured us that0\0is not equal to1, but rather

  • Ifa = 0, then0\0is Not-A-Number.

And, for a diagonal matrix of any order, this scalar case is applicable to each diagonal element.

Rank deficient matrices

IfAis arank deficientmatrix with rankr<n, thenA\Acannot possibly beI。The rank of the product of two matrices cannot be larger than rank of either matrix. SoA\Acannot outrankAitself and

  • IfAis rank deficient, thenA\Ais definitely notI

It just so happens that the most recent issue ofSIAM Reviewincludes the paper about matrix rank, "LU and CR Elimination", by my colleague Gil Strang and myself. The paper is available from either theSIAM web siteor Gil'sMIT web site。Another pointer is thisCleve's Corner

Lots of NaNs

Here are three examples whereA\AgeneratesNaN

A = 0 B = [1 0; 0 0] C = [1 2; 4 8] X = A\A Y = B\B Z = C\C
A = 0 B = 1 0 0 0 C = 1 2 4 8 X = NaN Y = 1 0 NaN NaN Z = NaN NaN NaN NaN

Magic squares

I always like to investigate any property of matrices by checking out magic squares.

small_fig warningoffnmax = 50; r = zeros(nmax,1); e = zeros(nmax,1);forn = 1:nmax A = magic(n); X = A\A; I = eye(n); r(n) = rank(A); e(n) = norm(X-I);end

Odd

MATLAB uses three different algorithms for computing magic squares, odd, singly even and doubly even. If the ordernis odd, thenA = magic(n)is nonsingular, the rank ofAisnand the elements of the computedA\Aare within roundoff error of the elements ofI。注意的程度r for the error plot is 3.0e-15.

n = 3:2:nmax; plots(n,r,60,e,[],"Odd")

Singly even

If the ordernis divisible by2, but not by4, thenmagic(n)is rank deficient. Its rank is about half its order. The error plot reflects the fact thatA\Ais notI

n = 2:4:nmax; plots(n,r,60,e,200,"Singly even")

Doubly even

If the ordernis divisible by4, thenmagic(n)is very rank deficient. The rank is always3。The error plots are all over the place. Orders8and40have errors that are larger than my plot scale. Orders16and32are missing entirely because computingA\Aencounters0\0resulting inNaN

n = 4:4:nmax; plots(n,r,12,e,750,"Doubly even")

Pseudoinverse

Ispinv(A)*bmore "robust" thanA\b?

You should not usepinvjust to create solutions to problems that do not have solutions. The pseudoinverse is intended to characterize the solution to a specific technical question: if a system of linear equations has many solutions, which is the shortest one? If you replaceA\bbypinv(A)*b, be sure that is what you want.

使用pinv代替反斜杠不做了ith rank deficiency. The difficulties are already present in mono-elemental matrices. The only rank deficient 1-by-1 matrix is0andpinv(0)=0。This is less in-your-face thanNaN, but there is no way you can makepinv(0)*0equal to1

When I redo the examples above, I get

A = 0 B = [1 0; 0 0] C = [1 2; 4 8] X = pinv(A)*A Y = pinv(B)*B Z = pinv(C)*C
A = 0 B = 1 0 0 0 C = 1 2 4 8 X = 0 Y = 1 0 0 0 Z = 0.2000 0.4000 0.4000 0.8000

The NaNs are gone, but is this really "more robust" than backslash? If you still think so, explain whereZcomes from.

My Bottom Line

This has been about square, dense, nonsymmetric matricesA。For such matrices:

  • A\Amay produceNaNfor rank-deficient matrices.
  • pinv(A)*AavoidsNaNby attempting to hide rank-deficient matrices.




Published with MATLAB® R2022a

|

コメント

コメントを残すには、ここをクリックして MathWorks アカウントにサインインするか新しい MathWorks アカウントを作成します。