1

In this question was found a solution to find a particular solution to a non-square linear system that has infinitely many solutions. This leads to another question:

How to find all the solutions for a non-square linear system with infinitely many solutions, with R? (see below for a possible description of the infinite set of solutions)

Example: the linear system 

x+y+z=1 
x-y-2z=2

is equivalent to A X = B with:

A=matrix(c(1,1,1,1,-1,-2),2,3,T)
B=matrix(c(1,2),2,1,T)

A
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1   -1   -2

B
     [,1]
[1,]    1
[2,]    2

We can describe the infinite set of solutions with:

x = 3/2 + (1/2) z
y = -1/2 + (-3/2) z
z in R

Thus, R could describe the set of solutions this way:

> solve2(A,B)
$principal             
[1] 1 2            # this means that x and y will be described

$free
[1] 3              # this means that the 3rd variable (i.e. z) is free in the set of real numbers

$P
[1] 1.5 -0.5

$Q
[1] 0.5 -1.5

This means that every solution can be created with:

z = 236782 # any value would be ok
solve2(A,B)$P + z * solve2(A,B)$Q     # this gives x and y

About the maths: there always exist such a decomposition, when the linear system has infinitely many solutions, this part is ok. The question is: is there something to do this in R?

Basj
  • 41,386
  • 99
  • 383
  • 673
  • This isn't clear to me. How do you know which is free? – kangaroo_cliff Oct 16 '17 at 09:29
  • Here the [rank](https://en.wikipedia.org/wiki/Rank_(linear_algebra)) of the linear system is 2, and there are 3 variables, thus 3-2=1 variable can be chosen as a free parameter. NB: to be more precise: "there exists 1 variable that can be chosen as a free paramter". Sometimes not *all of them* can be chosen, but at least one can. In general: it's possible to choose `number of variables - rank of the system` variables as free parameters (see linear algebra course). – Basj Oct 16 '17 at 09:31

2 Answers2

1

You can solve equations like thse using the generalized inverse of A.

library(MASS)
ginv(A) %*% B

# 1.2857143
# 0.1428571
#-0.4285714

A %*% ginv(A) %*% B

#    1
#    2

So, with help from @Bhas

gen_soln <- function(vec) {

  G <- ginv(A)
  W <- diag(3) - G %*% A
  (G %*% B + W %*% vec)
}

You can now find many solutions by providing a vector of length 3 to `gen_soln' function. For example, 

one_from_inf <- gen_soln(1:3)
one_from_inf
#[1,]  1.35714286
#[2,] -0.07142857
#[3,] -0.2857142

# Test the solution.

A %*% one_from_inf

#     [,1]
#[1,]    1
#[2,]    2

# Using random number generator 
A %*% gen_soln(rnorm(3))

#     [,1]
#[1,]    1
#[2,]    2
kangaroo_cliff
  • 6,067
  • 3
  • 29
  • 42
  • 1
    Read the section **Matrix Solution** on this page about solving [System_of_linear_equations](http://en.wikipedia.org/wiki/System_of_linear_equations) on Wikipedia. – Bhas Oct 16 '17 at 09:42
  • @Disco Superfly. Thanks for writing this explicitly. I was too lazy to spell it all out. I had the same method. A final remark: the solution `ginv(A) %*% B` is the minimum norm solution to the system of equations. – Bhas Oct 16 '17 at 12:55
1

The general solution to 

A*x = b

is

  x = x0 + z

where x0 is any solution and z is in the kernel of A

As pointed out above you can find a particular solution x0 by using the generalised inverse. You can also use the SVD to find a basis for the kernel of A:

A = U*S*V'

where U and V are orthogonal and S diagonal, with, say, the last k entries on the diagonal 0 (and the others non-zero).

If follows that the last k columns of V form a basis for the kernel of A, and if we call these z1,..zk then the solutions of the original equation are

x = x0 + c1*z1 + .. ck*zk

for any real c1..ck

dmuir
  • 4,211
  • 2
  • 14
  • 12
  • In my solution, I am actually using both the generalized inverse and the basis for kernel. I am not sure SVD is the best way to do it, for one it is never the last diagonal entries in S are 0 - they are quite small. – kangaroo_cliff Oct 17 '17 at 06:50
  • @Disco Yes, there are always inaccuracies in floating point calculations. You need to decide, both for computing the generalized inverse and the kernel basis which diagonal elements should be taken as zero. In the past I've used a (small!) multiple of the largest singular value as the measure of smallness, but have never been sure that this is the best thing to do – dmuir Oct 17 '17 at 08:25
  • I am not talking about the computational numerical errors. You can easily find non-square matrices with smallest singular value not close to zero; for example `svd(iris[11:13, -5])`. I guess one option would be to look at the % of variance explained. – kangaroo_cliff Oct 17 '17 at 12:11