Finding solution set of system of linear equations.

I have a simple system of linear equations to solve for X, aX=b:
a
[,1] [,2] [,3] [,4]
[1,]    1    2    1    1
[2,]    3    0    0    4
[3,]    1   -4   -2   -2
[4,]    0    0    0    0
b
[,1]
[1,]    0
[2,]    2
[3,]    2
[4,]    0

(This is ex Ch1, 2.2 of Artin, Algebra).
So, 3 eqs in 4 unknowns. One can easily use row-reductions to find a
homogeneous solution(b=0) of: 
X_1 = 0, X_2 = -c/2, X_3 = c,  X_4 = 0

and solutions of the above system are:
X_1 = 2/3, X_2 = -1/3-c/2,  X_3 = c, X_4 = 0.

So the Kernel is 1-D spanned by X_2 = -X_3 /2, (nulliity=1), rank is 3.

In R I use solve():
solve(a,b)
Error in solve.default(a, b) : 
  Lapack routine dgesv: system is exactly singular

and it gives the error that the system is exactly singular, since it seems
to be trying to invert a.
So my question is:
Can R only solve non-singular linear systems? If not, what routine should I
be using? If so, why? It seems that it would be simple and useful enough to
have a routine which, given a system as above, returns the null-space
(kernel) and the particular solution.

--
View this message in context: http://r.789695.n4.nabble.com/Finding-solution-set-of-system-of-linear-equations-tp3541490p3541490.html
Sent from the R help mailing list archive at Nabble.com.
solve() only works for nonsingular systems of equations.

Use a generalized inverse for singular systems:

 > A<- matrix(c(1,2,1,1, 3,0,0,4, 1,-4,-2,-2, 0,0,0,0), ncol=4, byrow=TRUE)
 > A
      [,1] [,2] [,3] [,4]
[1,]    1    2    1    1
[2,]    3    0    0    4
[3,]    1   -4   -2   -2
[4,]    0    0    0    0
 > b<- c(0,2,2,0)  #rhs
 > b
[1] 0 2 2 0
 >
 > require('MASS')
 > giA<- ginv(A) #M-P generalized inverse
 > giA
            [,1]          [,2]        [,3] [,4]
[1,]  0.6666667  1.431553e-16  0.33333333    0
[2,]  0.3333333 -1.000000e-01 -0.03333333    0
[3,]  0.1666667 -5.000000e-02 -0.01666667    0
[4,] -0.5000000  2.500000e-01 -0.25000000    0
 >
 > require('Matrix')
 > I<- as.matrix(Diagonal(4))  #order 4 identity matrix
 > I
      [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1
 >
 > giA%*%b   #particular solution
               [,1]
[1,]  6.666667e-01
[2,] -2.666667e-01
[3,] -1.333333e-01
[4,] -2.220446e-16
 > giA%*%A - I   #matrix for parametric homogeneous solution
               [,1] [,2] [,3]          [,4]
[1,]  0.000000e+00  0.0  0.0  5.551115e-16
[2,]  3.469447e-17 -0.2  0.4  4.024558e-16
[3,]  4.510281e-17  0.4 -0.8  2.706169e-16
[4,] -3.330669e-16  0.0  0.0 -7.771561e-16
I have a simple system of linear equations to solve for X, aX=b:
a
     [,1] [,2] [,3] [,4]
[1,]    1    2    1    1
[2,]    3    0    0    4
[3,]    1   -4   -2   -2
[4,]    0    0    0    0
b
     [,1]
[1,]    0
[2,]    2
[3,]    2
[4,]    0

(This is ex Ch1, 2.2 of Artin, Algebra).
So, 3 eqs in 4 unknowns. One can easily use row-reductions to find a
homogeneous solution(b=0) of:
X_1 = 0, X_2 = -c/2, X_3 = c,  X_4 = 0

and solutions of the above system are:
X_1 = 2/3, X_2 = -1/3-c/2,  X_3 = c, X_4 = 0.

So the Kernel is 1-D spanned by X_2 = -X_3 /2, (nulliity=1), rank is 3.

In R I use solve():
solve(a,b)
Error in solve.default(a, b) :
  Lapack routine dgesv: system is exactly singular

and it gives the error that the system is exactly singular, since it seems
to be trying to invert a.
So my question is:
Can R only solve non-singular linear systems? If not, what routine should I
be using? If so, why? It seems that it would be simple and useful enough to
have a routine which, given a system as above, returns the null-space
(kernel) and the particular solution.

--
View this message in context: 
http://r.789695.n4.nabble.com/Finding-solution-set-of-system-of-linear-equations-tp3541490p3541490.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
================================================================
Robert A. LaBudde, PhD, PAS, Dpl. ACAFS  e-mail: ral at lcfltd.com
Least Cost Formulations, Ltd.            URL: http://lcfltd.com/
824 Timberlake Drive                     Tel: 757-467-0954
Virginia Beach, VA 23464-3239            Fax: 757-467-2947

"Vere scire est per causas scire"
Thanks Robert. That all seems to work. I also found the MASS::Null() function
that gives the null space for the matrix(transpose) given as argument. I am
still trying to appreciate the math behind the Moore-Penrose inverse matrix.
If you have any suggestions for understanding how to use R to solve these
kinds of systems, please send me the pointer.

 I think the R documentation for solve() should state right up front that it
only solves non-singular systems, and it should point me to MASS::ginv() and
how to use it to solve this obvious kind of problem. Better, there could be
a generalized solve() which produces a particular solution, the null space,
image space (basis thereof). This other solution (using ginv() and giA%*%A -
I for the kernel) seems deeply embedded in a particular solution technique
(and should be available), but the generalized solve_lin_sys(), as I
suggested, seems generally quite useful..

Don

--
View this message in context: http://r.789695.n4.nabble.com/Finding-solution-set-of-system-of-linear-equations-tp3541490p3542930.html
Sent from the R help mailing list archive at Nabble.com.