Code covered by the BSD License

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

Highlights from Don't let that INV go past your eyes; to solve that system, FACTORIZE!

5.0
5.0 | 13 ratings Rate this file 22 Downloads (last 30 days) File Size: 610 KB File ID: #24119 Version: 1.8

Don't let that INV go past your eyes; to solve that system, FACTORIZE!

Tim Davis (view profile)

14 May 2009 (Updated )

A simple-to-use object-oriented method for solving linear systems and least-squares problems.

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

File Information
Description

What's the best way to solve a linear system? Backslash is simple, but the factorization it computes internally can't be reused to solve multiple systems (x=A\b then y=A\c). You might be tempted to use inv(A), as S=inv(A) ; x=S*b ; y=S*c, but's that's slow and inaccurate. The best solution is to use a matrix factorization (LU, CHOL, or QR) followed by forward/backsolve, but figuring out the correct method to use is not easy.
In textbooks, linear algebraic expressions often use the inverse. For example, the Schur complement is written as S=A-B*inv(D)*C. It is never meant to be used that way. S=A-B*(D\C) or S=A-(B/D)*C are better, but the syntax doesn't match the formulas in your textbook.

The solution is to use the FACTORIZE object. It overloads the backslash and mtimes operators (among others) so that solving two linear systems A*x=b and A*y=c can be done with this:

F=factorize(A) ;
x=F\b ;
y=F\c ;

An INVERSE method is provided which does not actually compute the inverse itself (unless it is explicitly requested). The statements above can be done equivalently with this:

S=inverse(A) ;
x=S*b ;
y=S*c ;

That looks like a use of INV(A), but what happens is that S is a factorized representation of the inverse of A, and multiplying S times another matrix is implemented with a forward/backsolve.

Wiith the INVERSE function, your Schur complement computation becomes S = A - B*inverse(D)*C which looks just like the equation in your book ... except that it computes a factorization and then does a forward/backsolve, WITHOUT computing the inverse at all.

An extensive demo is included, as well as a simple and easy-to-read version (FACTORIZE1) with fewer features, as an example of how to create objects in MATLAB.

If the system is over-determined, QR factorization is used to solve the least-squares problem. If the system is under-determined, the transpose of A is factorized via QR, and the minimum 2-norm solution is found. In MATLAB, x=A\b finds a "basic" solution to an under-determined system.

You can compute a subset of the entries of the inverse just by referencing them. For example, S=inverse(A), S(n,n) computes the (n,n) entry of the inverse, without computing all of them.

In the RARE case that you need the entire inverse (or pseduo-inverse of a full-rank rectangular matrix) you can force the object to become "double", with S = double (inverse (A)). This works for both full and sparse matrices (pinv(A) only works for full matrices).

And remember ... don't ever multiply a matrix by inv(A).

This code is published as "Algorithm 930: FACTORIZE: an object-oriented linear system solver for MATLAB", T. A. Davis, ACM Transactions on Mathematical Software, Vol 39, Issue 4, pp. 28:1 - 28:18, 2013.

Acknowledgements
Required Products MATLAB
MATLAB release MATLAB 7.12 (R2011a)
MATLAB Search Path
```/
/Factorize
/Factorize/Demo
/Factorize/Doc
/Factorize/Test```
Other requirements MATLAB 7.6 (R2008a) or later is required. SuiteSparse optionally required for sparse COD.
Tags for This File   Please login to tag files.
 Please login to add a comment or rating.
Comments and Ratings (18)
17 Mar 2016 MANIK BANSAL

MANIK BANSAL (view profile)

02 Feb 2016 Milind

Milind (view profile)

inv() is certainly slower than \ unless you have multiple right hand side vectors to solve for. But inv() is NOT inaccurate. The link elaborates further : http://arxiv.org/abs/1201.6035

>Several widely-used textbooks lead the >reader to believe that solving a >linear system of equations Ax = b by >multiplying the vector b by a computed >inverse inv(A) is inaccurate. >Virtually all other textbooks on >numerical analysis and numerical >linear algebra advise against using >computed inverses without stating >whether this is accurate or not. In >fact, under reasonable assumptions on >how the inverse is computed, x = >inv(A)*b is as accurate as the >solution computed by the best
>backward-stable solvers.

Comment only
23 Oct 2015 Eric Schols

Eric Schols (view profile)

Works as documented in the demo: like a charm!

One question: can this package also be used for eigendecomposition? For example:

A = rand(5);
[S,D] = eig(A);
norm(A-S*D*inverse(S))

But this still uses EIG, does this package contain a direct way to compute the eigendecomposition?

24 Apr 2015 Daniel

Daniel (view profile)

13 Nov 2014 Hamza

Hamza (view profile)

12 Oct 2011 Michael Völker

Michael Völker (view profile)

07 Sep 2011 Tim Davis

Tim Davis (view profile)

Comment from the author: I have addressed Ben's comment about uminus in this update (Sept 2011). You can now do -inverse(A)*b, or -2*inverse(A)*b, which you could not do in the previous version.

Comment only
22 Mar 2010 Tim Davis

Tim Davis (view profile)

-inverse(A)*b can also be done as -(inverse(A)*b), which is more natural than inverse(A)*(-b).

Comment only
12 Mar 2010 Tim Davis

Tim Davis (view profile)

Comment from the author: implementing the uminus function would be rather painful. I would need to return an object that is unchanged from the previous one, except with a flag stating that the result of all other functions must negate their result. This would be rather ugly and would needlessly complicate the code, just to add a rather minor feature.

Thanks for the suggestion ... I agree that uminus "ought" to work. It seems simple from the perspective of an end-user. But implementing it is far more trouble than it's worth, in my opinion.

Comment only
18 Jan 2010 Benjamin

Benjamin (view profile)

Very good, it is even faster when I solve a single linear system.

27 Sep 2009 Bruno Luong

Bruno Luong (view profile)

What mathworks are waiting to incorporate this package into Matlab?

05 Aug 2009 Ben

Ben (view profile)

All in all this package is a great idea and it looks like a first class implementation. As such, I'd like to help improve it by pointing out an issue. I am trying to solve Ax=-b, so I type in x=-inverse(A)*b, but I get the following error:

??? Undefined function or method 'uminus' for input arguments of type 'factorization_dense_lu'.

The work around is simple: x=inverse(A)*(-b). However, these are mathematically equivalent, so both should work.

Comment only
19 Jun 2009 Miao LI

Miao LI (view profile)

05 Jun 2009 Petter

Petter (view profile)

Nice with an update to modern MATLAB OOP.

27 May 2009 Petter

Petter (view profile)

This is an excellent piece of software. Makes great use of Matlab's OOP.

18 May 2009 Sebastien PARIS

Sebastien PARIS (view profile)

Great !!!

15 May 2009 John D'Errico

John D'Errico (view profile)

Absolutely splendid! Many thanks to Tim for writing this. A great title too.

15 May 2009 us

us (view profile)

a truly professional and INValuable package written with a twinkle in the (author's) eye and a must look-at for people dealing with linear systems...
one pedestrian thought: why not create a stand-alone object given ML's new and powerful OOP-strength...
us

Updates
18 May 2009 1.1

New overload for plus and minus (cholupdate). Improved performance. More extensive demo.

19 May 2009 1.2

Bug fix so that inverse(A)\b does the same thing as A*b (a peculiar case, of course...)

20 May 2009 1.3

Added "(SetAccess = protected)" to class properties.

27 May 2009 1.4

Added link to LINFACTOR file.

04 Jun 2009 1.5

Substantial update. The factorization class is now a superclass, with 8 specific subclasses beneath it (one for each kind of matrix factorization). Minor changes to the demo. No change to the user interface.

06 Sep 2011 1.7

Major update. Added SVD and COD (complete orthogonal decomposition) and many methods, such as the uminus requested in a comment below.

20 Nov 2014 1.8

Converted to a toolbox. Updated to the version that appears as Algorithm 930 in the Collected Algorithms of the ACM.

Contact us