Don't Invert the Matrix

Original Post - 18 Nov 2022 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

A common issue with linear algebra textbooks is the depiction of x=A^-1b as the solution to the linear system of equations Ax=b.

Find the inverse of the matrix A, then multiply that inverse with the right-hand side vector, b.

Theoretically correct? Yes. Practical? No. Numerical linear algebraicists? They scream.

The practical approach to find x is direct solution by Gauss elimination (LU decomposition). In MATLAB speak, this operation is x=A \ b, a notation that I will adopt in the equation environments of my lecture notes.

Not only is inversion and multiplication more expensive than direct solution for x, it is also less accurate. Google something like “reasons to not invert a matrix” and you’ll find plenty of explanations, such as this one.

There’s a few element and material models in OpenSees that use the “by the book” inversion and multiplication approach. Usually, the offending operations appear in a local Newton loop where the Matrix::Invert() method is called in preparation for multiplication on the next line.

Matrix A(4,4); // Jacobian
Matrix Ainv(4,4); // Inverse of Jacobian
Vector x(4); // Unknowns
Vector b(4); // Right-hand side

// Form initial b

// Solve using Newton's method
while (b.Norm() > tol) {

    //
    // Form A based on current x
    //

    A.Invert(Ainv); // DON'T DO THIS
    dx = Ainv*b;
    x -= dx;

    //
    // Update b based on current x
    //
}

The use of the overloaded * operator to get dx is an offense that I’m willing to overlook.

What you should do instead is call the Matrix::Solve() method to get dx directly.

    A.Solve(b,dx);
    x -= dx;

You may say the matrices at the element and material level are pretty small, so it’s not a big deal. But all the unnecessary computations add up.

Also, don’t get me wrong, there are some instances where computing and storing the inverse is necessary within some element and material models in OpenSees, e.g., in the mixedBeamColumn element. In this formulation, a matrix inverse is used a few times across multiple operations. Ideally, we’d have the Matrix class store the factors from direct solution so that one could call Matrix::Solve() multiple times, performing backward substitution after the first call.

But that’s a PR for another day. And Frank doesn’t like it when people mess with the Matrix class.

OpenSees Cloud

OpenSees AMI

Don't Invert the Matrix