Purpose of inv
I have received comments from several savvy customers suggesting that we remove the function inv from MATLAB. The reasons are because people who don't know enough numerical analysis will use it blindly and incorrectly and often get a poor quality answer unnecessarily.
Contents
Reasons to Keep inv
There are several reasons why we do keep the function inv. These reasons are:
- for pedagogical and teaching use.
- historical (we'd create an incompatibility).
- meeting expectations. This is somewhat similar to the factorial, function; users expected it to be there even though they could use gamma instead).
- for certain computations. An example is when you need the covariance of least squares estimates. It's hard to accomplish that without inv.
Reasons to Not Use or Remove inv
There are good reasons to not use inv however. The main one is
- numerical considerations.
Here is the second paragraph of the description for inv in the documentation:
In practice, it is seldom necessary to form the explicit inverse of a matrix. A frequent misuse of inv arises when solving the system of linear equations.
One way to solve this is with x = inv(A)*b. A better way, from both an execution time and numerical accuracy standpoint, is to use the matrix division operator x = A\b. This produces the solution using Gaussian elimination, without forming the inverse. See \ and / for further information.
According to Cleve, the best way to compute inv(X'*X) is to avoid computing X'*X, and use the economy-sized QR factorization of X instead.
[Q,R] = qr(X,0)
S = inv(R)
S*S'
To see an extreme example, try
d = 1.e-10 X = [1 1; d 0; 0 d]
d =
1.0000e-010
X =
1.0000 1.0000
0.0000 0
0 0.0000
inv(X'*X)
Warning: Matrix is singular to working precision. ans = Inf Inf Inf Inf
[Q,R] = qr(X,0) S = inv(R) S*S'
Q =
-1.0000 0.0000
-0.0000 -0.7071
0 0.7071
R =
-1.0000 -1.0000
0 0.0000
S =
1.0e+009 *
-0.0000 -7.0711
0 7.0711
ans =
1.0e+019 *
5.0000 -5.0000
-5.0000 5.0000
The computed X' * X is exactly singular and does not have an inverse, but R does have an inverse. We're still using the inv function, but on a nonsingular matrix.
Your Thoughts?
So, we have the function inv that is ripe for misuse. What's the best way we can steer users to better solutions? Please post your thoughts here.
