0

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?

Kristoffer Lindvall
  • 121
  • 2
  • 11

3 Answers3

2

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]
divenex
  • 15,176
  • 9
  • 55
  • 55
andersource
  • 799
  • 3
  • 9
  • One should compute `x_qr = np.linalg.lstsq(A, b)[0]` directly, without first computing the QR decomposition of `A` which is redundant and is already done by `lstsq`. – divenex Dec 11 '19 at 14:41
1

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).

Tunneller
  • 381
  • 2
  • 13
0

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)
yatu
  • 86,083
  • 12
  • 84
  • 139