CODEWITHSUNDEEP
pythonnumpy
Jan Erst
Jan Erst

Reputation: 89

numpy linalg.solve, not a square matrix

So currently I'm working with code looking like:

Q,R = np.linalg.qr(matrix)
Qb = np.dot(Q.T, new_mu[b][n])
x_qr = np.linalg.solve(R, Qb)
mu.append(x_qr)

The code works fine as long as my matrix is square, but as soon as it's not, the system is not solvable and I got errors. If I've understood it right I can't use linalg.solve on non-full rank matrices, but is there a way for me to get across this obstacle without using a lstsquare solution?

Upvotes: 3

Views: 4757

Answers (2)

JeeyCi
JeeyCi

Reputation: 579

you should try to incorporate rank-revealing QR decomposition with column-pivoting (e.g. like here) to below-mentioned solve or Normal Equation solution, lstsq as I remember, probably, uses SVD-decomposition to get rank of linear system (that in general is considered to result in more numerically stable solutions) - remember rank-nullity theorem

p.s.

from LS Normal Equation intercept, *coef = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y) we can use equation (X.T.dot(X)).dot(β) = (X.T).dot(y)

from  scipy import linalg
import numpy as np
from sklearn.datasets import make_regression

# Generate a regression problem
X, y = make_regression(
    n_samples=50,
    n_features=2,
    n_informative=2,
    noise = 10,
    random_state=25
    )

X_b = np.column_stack((np.ones(len(y)), X)) # for intercept

########################## solve for non-square
XtX = X_b.T.dot(X_b)
XtY = X_b.T.dot(y)
x1 = np.linalg.solve(XtX, XtY);
print(x1)

########################## least-squares - compute vector x such that 2-norm is minimized
x2, sse, rank, sva = linalg.lstsq(X_b, y)
print(x2)

########### check
# euclidean distance - here shows x_solve ~ x_lstsq
print("distance between vectors: ", np.linalg.norm(x1-x2,2))
print(np.allclose(x1, x2))

## distance between vectors:   - as is tiny difference
## True

######################### check with ML
from sklearn.linear_model import LinearRegression

lr = LinearRegression()
lr.fit(X, y)

print(f"Intercept: {lr.intercept_}\n\
Coefficients: {lr.coef_}")

Upvotes: 0

Kevin Ji
Kevin Ji

Reputation: 10489

No, this is not possible, as specified in the np.linalg.solve docs.

The issue is that given Ax = b, if A is not square, then your equation is either over-determined or under-determined, assuming that all rows in A are linearly independent. This means that there does not exist a single x that solves this equation.

Intuitively, the idea is that if you have n (length of x) variables that you are trying to solve for, then you need exactly n equations to find a unique solution for x, assuming that these equations are not "redundant". In this case, "redundant" means linearly dependent: one equation is equal to the linear combination of one or more of the other equations.

In this scenario, one possibly useful thing to do is to find the x that minimizes norm(b - Ax)^2 (i.e. linear least squares solution):

x, _, _, _ = np.linalg.lsq(A, b)

Upvotes: 5

Related Questions