From c3ce78725469fb9d75d7c1d74b33bf1d65eea9cc Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sun, 4 Oct 2020 08:27:37 -0400 Subject: [PATCH] mjo/ldlt.py: rule #1 is never compute the inverse of a matrix. --- mjo/ldlt.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/mjo/ldlt.py b/mjo/ldlt.py index 6db747f..9a6070d 100644 --- a/mjo/ldlt.py +++ b/mjo/ldlt.py @@ -537,11 +537,17 @@ def block_ldlt(A): C = A[k+2:n,k:k+2] B = A[k+2:,k+2:] - # TODO: don't invert, there are better ways to get the C*E^(-1) - # that we need. - E_inverse = E.inverse() + # We don't actually need the inverse of E, what we really need + # is C*E.inverse(), and that can be found by setting + # + # C*E.inverse() == X <====> XE == C. + # + # The latter can be found much more easily by solving a system. + # Note: I do not actually know that sage solves the system more + # intelligently, but this is still The Right Thing To Do. + CE_inverse = E.solve_left(C) - schur_complement = B - (C*E_inverse*C.transpose()) + schur_complement = B - (CE_inverse*C.transpose()) # Compute the Schur complement that we'll work on during # the following iteration, and store it back in the lower- @@ -562,7 +568,7 @@ def block_ldlt(A): # the lower-left-hand corner of "A", being sure to set # the lower-left entries from the upper-right ones to # avoid collisions. - A[k+i+2,k+j] = (C*E_inverse)[i,j] + A[k+i+2,k+j] = CE_inverse[i,j] k += 2 -- 2.44.2