From: Michael Orlitzky <michael@orlitzky.com>
Date: Sun, 4 Oct 2020 12:27:37 +0000 (-0400)
Subject: mjo/ldlt.py: rule #1 is never compute the inverse of a matrix.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=c3ce78725469fb9d75d7c1d74b33bf1d65eea9cc;p=sage.d.git

mjo/ldlt.py: rule #1 is never compute the inverse of a matrix.
---

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