Search code examples
pythonnumpyleast-squaresqr-decomposition

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?

Solution

  • 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]