Solve linear equations on a Gram matrix with numpy

0
votes

Consider a case where, given an MxM matrix A and a vector b, I want to solve something of the form inv(A @ A.T) @ b (where I know A is invertible). As far as I know, it is always faster to use solve_* rather than inv. There are also variants for more efficient solving for PSD matrices (which A @ A.T must be), using Cholesky factorization.

My question - since I'm constructing the matrix A @ A.T just to immediately throw it away - is there a more specialized procedure for solving linear equations with the gram matrix of A without having to construct it?

1
pythonnumpyscipylinear-algebralapack
This sounds like it would be a good question to ask over in scicomp.stackexchange.com - Warren Weckesser

1 Answers

2
votes

You can compute the factorization of A and then use that to solve your system.

Assume we want to solve

A A^T x = b

for x.

Compute the factorization of A=LU. Then solve Ay=b for y. Then solve A^T x = y for x.

This way you dont have to compute the matrix A^T A.

Note that if one has a factorization of A=LU then one can solve Ax=b as well as A^T x=b efficiently for x. This is because A^T=U^T L^T which is again a factorization of a lower times an upper triangular matrix.