Solving Non-Square Linear System with R

Solving non-square linear system with R

If the matrix A has more rows than columns, then you should use least squares fit.

If the matrix A has fewer rows than columns, then you should perform singular value decomposition. Each algorithm does the best it can to give you a solution by using assumptions.

Here's a link that shows how to use SVD as a solver:

http://www.ecse.rpi.edu/~qji/CV/svd_review.pdf

Let's apply it to your problem and see if it works:

Your input matrix A and known RHS vector B:

> A=matrix(c(0,1,-2,3,5,-3,1,-2,5,-2,-1,1),3,4,T)
> B=matrix(c(-17,28,11),3,1,T)
> A
     [,1] [,2] [,3] [,4]
[1,]    0    1   -2    3
[2,]    5   -3    1   -2
[3,]    5   -2   -1    1
> B
     [,1]
[1,]  -17
[2,]   28
[3,]   11

Let's decompose your A matrix:

> asvd = svd(A)
> asvd
$d
[1] 8.007081e+00 4.459446e+00 4.022656e-16

$u
           [,1]       [,2]       [,3]
[1,] -0.1295469 -0.8061540  0.5773503
[2,]  0.7629233  0.2908861  0.5773503
[3,]  0.6333764 -0.5152679 -0.5773503

$v
            [,1]       [,2]       [,3]
[1,]  0.87191556 -0.2515803 -0.1764323
[2,] -0.46022634 -0.1453716 -0.4694190
[3,]  0.04853711  0.5423235  0.6394484
[4,] -0.15999723 -0.7883272  0.5827720

> adiag = diag(1/asvd$d)
> adiag
          [,1]      [,2]        [,3]
[1,] 0.1248895 0.0000000 0.00000e+00
[2,] 0.0000000 0.2242431 0.00000e+00
[3,] 0.0000000 0.0000000 2.48592e+15

Here's the key: the third eigenvalue in d is very small; conversely, the diagonal element in adiag is very large. Before solving, set it equal to zero:

> adiag[3,3] = 0
> adiag
          [,1]      [,2] [,3]
[1,] 0.1248895 0.0000000    0
[2,] 0.0000000 0.2242431    0
[3,] 0.0000000 0.0000000    0

Now let's compute the solution (see slide 16 in the link I gave you above):

> solution = asvd$v %*% adiag %*% t(asvd$u) %*% B
> solution
          [,1]
[1,]  2.411765
[2,] -2.282353
[3,]  2.152941
[4,] -3.470588

Now that we have a solution, let's substitute it back to see if it gives us the same B:

> check = A %*% solution
> check
     [,1]
[1,]  -17
[2,]   28
[3,]   11

That's the B side you started with, so I think we're good.

Here's another nice SVD discussion from AMS:

http://www.ams.org/samplings/feature-column/fcarc-svd

Solve non square linear system with Math.net

You can not use the Matrix.Solve function with a non-square matrix because there is no inverse and no unique solutions for a rectangular matrix. Google "inverse of rectangular matrix" for explanations galore. You can use pseudoinverse however, as shown below.

        var mBuilder = Matrix<double>.Build;
        var A = mBuilder.DenseOfArray(new double[,]
      {
          { 3, 2, 1, 5, -1, 0, 0 },
          { 2, 1, 1, 2, 0, -1, 0 },
          { 5, 1, 3, 4, 0, 0, -1 }
      });

        Matrix<double> b = Matrix<double>.Build.Dense(3, 1);
        b[0, 0] = -3.0;
        b[1, 0] = -5.0;
        b[2, 0] = -2.0;

        var p = A.PseudoInverse();
        var x = p * b;

        // verify
        var o = A * x;

Solving non-square system of linear equations in Javascript

After a lot of trial and error I found out that the solve function in ml-matrix allows for more than N equations with N variables.

I also tried the approach with pseudo-inverse, but that also required the matrix to be square.

Solving a system of linear equations in a non-square matrix

Whether or not your matrix is square is not what determines the solution space. It is the rank of the matrix compared to the number of columns that determines that (see the rank-nullity theorem). In general you can have zero, one or an infinite number of solutions to a linear system of equations, depending on its rank and nullity relationship.

To answer your question, however, you can use Gaussian elimination to find the rank of the matrix and, if this indicates that solutions exist, find a particular solution x0 and the nullspace Null(A) of the matrix. Then, you can describe all your solutions as x = x0 + xn, where xn represents any element of Null(A). For example, if a matrix is full rank its nullspace will be empty and the linear system will have at most one solution. If its rank is also equal to the number of rows, then you have one unique solution. If the nullspace is of dimension one, then your solution will be a line that passes through x0, any point on that line satisfying the linear equations.

Is there a way to generate random solutions to non-square linear equations, preferably in python?

Thanks for the answers, which both helped me and pointed me in the right direction.

I now have an easy step-by-step solution to my problem for arbitrary k<n.

1. Find one solution to all equations given. This can be done by using

 solution_vec = numpy.linalg.lstsq(A,b)

this gives a solution as seen in ukemis answer. In my example above, the Matrix A is equal to the coefficients of the equations on the left side, b represents the vector on the right side.

2. Determine the null space of your matrix A.

These are all vectors v such that the skalar product v*A_i = 0 for every(!) row A_i of A. The following function, found in this thread can be used to get representatives of the null space of A:

def nullSpaceOfMatrix(A, eps=1e-15):
    u, s, vh = scipy.linalg.svd(A)
    null_mask = (s <= eps)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)  

3. Generate as many (N) "random" linear combinations (meaning with random coefficients) of solution_vec and resulting vectors of the nullspace of the matrix as you want! This works because the scalar product is additive and nullspace vectors have a scalar product of 0 to the vectors of the equations. Those linear combinations always must contain solution_vec, as in:

linear_combination = solution_vec + a*null_spacevec_1 + b*nullspacevec_2...

where a and b can be randomly chosen.

Solving a simply non linear system

R can do this quite easily

library(nleqslv)
f <- function(x) {
  a<-x[1]
  b<-x[2]
  c(a*b-m1,a*b^2+(a^2)*(b^2)-m2)
}
nleqslv(c(1,30), f)

Output should look like:

$`x`
[1]  1.286486 30.595840

$fvec
[1] -9.663381e-13 -1.396074e-10

$termcd
[1] 1

$message
[1] "Function criterion near zero"

$scalex
[1] 1 1

$nfcnt
[1] 11

$njcnt
[1] 2

$iter
[1] 10

You can make things more robust by providing gradients. Of course R can also estimate parameters for the gamma distribution directly (e.g. fitdistr from the MASS package).

---

Related Topics