9

I'm focusing on the special case where A is a n x d matrix (where k < d) representing an orthogonal basis for a subspace of R^d and b is known to be inside the subspace. I thought of using the tools provided with numpy, however they only work with square matrices. I had the approach of filling in the matrix with some linearly independent vectors to "square" it and then solve, but I could not figure out how to choose those vectors so that they will be linearly independent of the basis vectors, plus I think it's not the only approach and I'm missing something that can make this easier. is there indeed a simpler approach than the one I mentioned? if not, how indeed do I choose those vectors that would complete Ainto a square matrix?

Ron Tubman
  • 99
  • 1
  • 1
  • 4

4 Answers4

13

As you mentioned, np.linalg.solve needs a full-rank square matrix.

For all the other linear cases, if you are interested in min||Ax-b||^2. (you probably are) you can use np.linalg.lstsq.

Solves the equation a x = b by computing a vector x that minimizes the Euclidean 2-norm || b - a x ||^2.

The equation may be under-, well-, or over- determined (i.e., the number of linearly independent rows of a can be less than, equal to, or greater than its number of linearly independent columns). If a is square and of full rank, then x (but for round-off error) is the “exact” solution of the equation.

(bold annotation by me)

This is also mentioned in the docs of the original np.linalg.solve:

a must be square and of full-rank, i.e., all rows (or, equivalently, columns) must be linearly independent; if either is not true, use lstsq for the least-squares best “solution” of the system/equation.

Community
  • 1
  • 1
sascha
  • 32,238
  • 6
  • 68
  • 110
  • 1
    What if I want the solution to satisfy a given inequality such as all coefficients should be positive? – gota Mar 09 '18 at 12:54
  • 1
    @gota Then you need constraint-optimization methods, which are more complex. As there is no closed-form solution anymore all of these are iterative. scipy.optimize has some candidates, depending on the task. In the case of var-lower-bounds: either optimize.nnls (dense, small) or a customized optimize.minimize(method='L-BFGS-B') with bounds. – sascha Mar 09 '18 at 15:38
3

If you have fewer equations than unknowns (assuming you meant to type n < d), you are not expecting a unique solution. You can obtain the solutions using a singular value decomposition.

  • Start by finding the SVD, consisting of three matrices U S V^T using numpy.linalg.svd. Then you can get the pseudoinverse of A from V [ diag(1/s_j) ] U^T, which gives you one solution.
  • Your final solution will be the one solution you found, plus linear combinations of nullspace basis vectors of A.
  • To find the nullspace basis vectors of A, extract columns j from the matrix V that correspond to singular values s_j from the matrix S that are zero (or, below some "small" threshold).

You can implement this last bit pretty easily in Python using for loops & if statements - the heavy lifting is the decomposition itself. Press et al's excellent Numerical Recipes covers linear algebra in Chapter 2 (freely available version of Numerical Recipes in C here and here). They have a great presentation of SVD that explains both the theory of SVD and how that translates into algorithms (mainly focusing on how to use the result of an SVD). They provide an SVD object that has more functionality than numpy.linalg.svd, but as mentioned, the core functionality is the actual decomposition, and doing things like getting nullspace basis vectors is just dressing around for loops and if statements operating on U, S, and V.

charlesreid1
  • 4,360
  • 4
  • 30
  • 52
  • Can you please clarify how does getting the pseudoinverse give me one solution? and in addition, I hurriedly forgot to mention that the rows of A are known to be orthogonal. does that make the entire thing easier? – Ron Tubman Sep 23 '17 at 10:41
  • In addition, if I do have an infinite amount of solutions, how can I find those (if any) that satisfy a given constraint – Ron Tubman Sep 23 '17 at 11:04
  • You can find a solution using the pseudoinverse of A the same way you would use the inverse of A: X = pseudoinverse(A) * B. If A is already row-orthogonal, the algorithm would be easier, but you'd have to modify the SVD algorithm itself. SVD algorithms can be pretty complicated, so probably best just to stick with using the Numpy implementation and sacrifice any potential savings. – charlesreid1 Sep 23 '17 at 19:29
  • Regarding your second question - SVD can tell you if your matrix is singular (if all the singular values are zero, your matrix is singular & you have an infinite number of solutions). But it definitely can't solve the problem of solving a linear system under constraints - you probably want something in the [scipy.optimize](https://docs.scipy.org/doc/scipy-0.15.1/reference/optimize.html) package. – charlesreid1 Sep 23 '17 at 19:38
  • Can you post a small example so we could make this a bit more concrete? – charlesreid1 Sep 23 '17 at 19:39
  • let us say I'm working in 4-dimensional space, and I have these 2 vectors: v_1 = [2,6,3,7], v_2 = [5, 8, 1, 4]. these are not necessarily orthogonal, so I generated (in the past) an orthogonal basis for them using, say, 'numpy.linalg.qr'. those vectors span some subspace, and I have now a new vector, let's say, u = [9, 20, 7, 18], that is on that space (Let's say I know it because I could compute the distance to that space using projection on the Orthonormal basis I computed), and I want it's coordinates in that space, that is the coefficients a_2, a_2 such as u = a_1 v_1 + a_2 v_2. – Ron Tubman Sep 24 '17 at 15:00
0

Can use:

  1. np.linalg.lstsq(x, y)
  2. np.linalg.pinv(X) @ y
  3. LinearRegression().fit(X, y)

Where 1 and 2 are from numpy and 3 is from sklearn.linear_model.

As a side note you will need to concatenate ones(use np.ones_like) in both 1 and 2 to represent the bias from the equation y = ax + b

Diogo Santiago
  • 229
  • 2
  • 9
0

In addition to the Numpy resources mentioned above, I would like to add the Solveset module from SymPy, and in particular linsolve. As stated in their documentation:

Solve system of N linear equations with M variables; both underdetermined and overdetermined systems are supported. The possible number of solutions is zero, one or infinite. Zero solutions throws a ValueError, whereas infinite solutions are represented parametrically in terms of the given symbols. For unique solution a FiniteSet of ordered tuples is returned.

They go on to provide examples and describe the supported formats for inputting Ax = b. Personally, I have had a lot of success incorporating this module into my Python projects, so I wanted to add it to the conversation.

octocat
  • 1
  • 1