Geek Answers Handbook

Getting a reversible square matrix from a non-square matrix of full rank in numpy or matlab

Suppose you have a full-rank matrix NxM A , where M>N If we denote the columns on C_i (with dimensions Nx1 ), we can write the matrix as

 A = [C_1, C_2, ..., C_M]

How can you get the first linearly independent columns of the original matrix A , so that you can build a new matrix NxN B , which is an invertible matrix with a nonzero determinant.

 B = [C_i1, C_i2, ..., C_iN]

How can you find indices {i1, i2, ..., iN} either matlab or python numpy? Is it possible to do this by expanding on singular values? Code snippets will be very welcome.

EDIT: To make this more specific, consider the following python code

 from numpy import * from numpy.linalg.linalg import det M = [[3, 0, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 0, 1], [0, 2, 0, 0, 0]] M = array(M) I = [0,1,2,4] assert(abs(det(M[:,I])) > 1e-8)

Therefore, given the matrix M, we need to find the indices of the set N linearly independent column vectors.

+4
source share
2 answers

Easy, peasy in MATLAB. Use a QR, in particular a hinged QR.

 M = [3 0 0 0 0; 0 0 1 0 0; 0 0 0 0 1; 0 2 0 0 0] [Q,R,E] = qr(M) Q = 1 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 R = 3 0 0 0 0 0 2 0 0 0 0 0 1 0 0 0 0 0 1 0 E = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0

The first 4 columns of E indicate the columns of M to be used, that is, the columns [1,2,3,5]. If you need M columns, just form the product M * E.

 M*E ans = 3 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 2 0 0 0

By the way, using det to determine if a matrix is ​​singular is absolutely, positively, definitely the worst way you can do this.

Use rank instead.

In fact, you should NEVER use det in MATLAB unless you understand why this is so bad, and you decide to use it, despite this fact.

+6
source

My first thought is to try every possible combination of N from the columns of M. This can be done like this (in Python):

 import itertools import numpy.linalg # 'singular' returns whether a matrix is singular. # You could use something more efficient than the determinant # (I'm not sure what options there are in NumPy) singular = lambda m: numpy.linalg.det(m) == 0 def independent_square(A): N,M = A.shape for colset in itertools.combinations(xrange(M), N): B = A[:,colset] if not singular(B): return B

If you need column indices instead of the resulting square matrix, just replace return B with return colset . Or you can get both using return colset,B

I do not know how SVD will help here. In fact, I can’t come up with any purely mathematical operation that would convert A to B (or even any one that could determine the column selection matrix MxN Q, such that B = AQ), with the exception of trial and error. But if you want to find out if one exists, math.stackexchange.com will be a good place to ask.

If all you need is a way to do it computationally, the above code should be sufficient.

+1
source

More articles:


All Articles