CODEWITHSUNDEEP
pythonnumpyleast-squaresqr-decomposition
Kristoffer Lindvall
Kristoffer Lindvall

Reputation: 141

Solve overdetermined system with QR decomposition in Python

I'm trying to solve an overdetermined system with QR decomposition and linalg.solve but the error I get is

LinAlgError: Last 2 dimensions of the array must be square.

This happens when the R array is not square, right? The code looks like this

import numpy as np
import math as ma

A = np.random.rand(2,3)
b = np.random.rand(2,1) 
Q, R = np.linalg.qr(A)
Qb = np.matmul(Q.T,b)
x_qr = np.linalg.solve(R,Qb)

Is there a way to write this in a more efficient way for arbitrary A dimensions? If not, how do I make this code snippet work?

Upvotes: 1

Views: 4727

Answers (4)

Royi
Royi

Reputation: 4963

For updated solution for the case solving multiple right hand sides for the same non square matrix:

import numpy as np
import scipy as sp

numRows = 500
numCols = 100
numIn   = 1000 #<! Number of inputs

mA = np.random.randn(numRows, numCols)
mB = np.random.randn(numRows, numIn)

mX = np.zeros(shape = (numCols, numIn))

mQ, mR = sp.linalg.qr(mA, mode = 'economic')

for ii in range(numIn):
    mX[:, ii] = sp.linalg.solve_triangular(mR, mQ.T @ mB[:, ii], check_finite = False)

This should be faster than calling higher level functions like solve() or lstsq().

Upvotes: 0

Tunneller
Tunneller

Reputation: 589

You need to call QR with the flag mode='reduced'. The default Q R matrices are returned as M x M and M x N, so if M is greater than N then your matrix R will be nonsquare. If you choose reduced (economic) mode your matrices will be M x N and N x N, in which case the solve routine will work fine.

However, you also have equations/unknowns backwards for an overdetermined system. Your code snippet should be

import numpy as np 

A = np.random.rand(3,2)
b = np.random.rand(3,1) 
Q, R = np.linalg.qr(A, mode='reduced')
#print(Q.shape, R.shape)
Qb = np.matmul(Q.T,b)
x_qr = np.linalg.solve(R,Qb)

As noted by other contributors, you could also call lstsq directly, but sometimes it is more convenient to have Q and R directly (e.g. if you are also planning on computing projection matrix).

Upvotes: 2

andersource
andersource

Reputation: 829

The reason is indeed that the matrix R is not square, probably because the system is overdetermined. You can try np.linalg.lstsq instead, finding the solution which minimizes the squared error (which should yield the exact solution if it exists).

import numpy as np

A = np.random.rand(2, 3)
b = np.random.rand(2, 1) 
x_qr = np.linalg.lstsq(A, b)[0]

Upvotes: 4

yatu
yatu

Reputation: 88236

As shown in the documentation of numpy.linalg.solve:

Computes the “exact” solution, x, of the well-determined, i.e., full rank, linear matrix equation ax = b.

Your system of equations is underdetermined not overdetermined. Notice that you have 3 variables in it and 2 equations, thus fewer equations than unknowns.

Also notice how it also mentions that in numpy.linalg.solve(a,b), a must be an MxM matrix. The reason behind this is that solving the system of equations Ax=b involves computing the inverse of A, and only square matrices are invertible.

In these cases a common approach is to take the Moore-Penrose pseudoinverse, which will compute a best fit (least squares) solution of the system. So instead of trying to solve for the exact solution use numpy.linalg.lstsq:

x_qr = np.linalg.lstsq(R,Qb)

Upvotes: 1

Related Questions