CODEWITHSUNDEEP
python-3.xnumpyscipylinear-algebra
Mostafa Orooji
Mostafa Orooji

Reputation: 81

numpy.linalg.solve with multiple dependent right-hand side

I'm trying to solve an equation system like Ax=b with full rank matrix A and right-hand side b including two columns for two different systems with the same left side. (Ax=b1, Ax=b2, b=concatenate((b1, b2), axis=1)). The point is that we have the same A for two systems. Therefore we should be able to use info from the first system for the second i.e the inverse of A.

numpy.linalg.solve easily solves this system if the columns of b are independent that is not my case. In my case, the second column of b depends on the solution of the first system.

I have tried to factorize matrix A and use this factorization to solve two systems. But, it is not efficient at all. Considering A is not symmetric I have used RQ and LU decomposition. 

from scipy.linalg import lu
import numpy as np
from scipy.linalg import solve_triangular
def qr_solver(a, b):
   q, r = np.linalg.qr(a)
   z = solve_triangular(r, b, lower=False)
   ans = np.dot(q.T, z)
   return ans
def plu_solver(a, b):
   per, l, u = lu(kkt_matrix)
   z = np.dot(per.T , b)
   x = solve_triangular(l, z, lower=True)
   ans = solve_triangular(u, x)
   return ans

Upvotes: 0

Views: 1485

Answers (2)

Mostafa Orooji
Mostafa Orooji

Reputation: 81

I just tried to just program what 'Talonmies' stated:

    from scipy.linalg import lu_factor, lu_solve
    from scipy.linalg import cho_factor
    from scipy.linalg import solve
    import numpy as np
    import time
    # Solving Ax = b1, Ay = f(x) with same A
    A = np.random.normal(1, 10, (2000, 2000))
    b1 = np.random.normal(1, 10, (2000, 1))
    b2 = np.random.normal(1, 10, (2000, 1))
    start = time.time()
    lu, pivot = lu_factor(A)
    x = lu_solve((lu, pivot), b1)
    y = lu_solve((lu, pivot), b2)
    end = time.time()
    time_factorization = start-end
    start = time.time()
    x = np.linalg.solve(A, b1)
    y = np.linalg.solve(A, b1)
   end = time.time()
   time_solve = start-end
   print('time_factorization:', time_factorization, 'time_solve:', time_solve)

and here the results: time_factorization: 0.6500372886657715 time_solve: 0.9420537948608398

It seems that this kind of factorization(getting lu and pivots) is really efficient.

Upvotes: 0

talonmies
talonmies

Reputation: 72349

scipy exposes lu_factor and lu_solve for this sort of problem:

from scipy.linalg import lu_factor, lu_solve

# Solving Ax = b1, Ay = f(x) with same A

lu, pivot = lu_factor(A)
x = lu_solve((lu, pivot), b1)
b2 = f(x)
y = lu_solve((lu, pivot), b2)

So if the RHS vectors are not linearly independent (implicit Runge-Kutta schemes are a good example), you can factorize the LHS once, and re-use it to solve as often as required.

Upvotes: 2

Related Questions