> On 20 Apr 2016, at 13:22, A A via R-help <r-help@r-project.org> wrote:
> 
> 
> 
> 
> I have a situation in R where I would like to find any x (if one exists) that 
> solves the linear system of equations Ax = b, where A is square, sparse, and 
> singular, and b is a vector. Here is some code that mimics my issue with a 
> relatively simple A and b, along with three other methods of solving this 
> system that I found online, two of which give me an error and one of which 
> succeeds on the simplified problem, but fails on my data set(attached). Is 
> there a solver in R that I can use in order to get x without any errors given 
> the structure of A? Thanks for your time.
> #CODE STARTS HEREA = 
> as(matrix(c(1.5,-1.5,0,-1.5,2.5,-1,0,-1,1),nrow=3,ncol=3),"sparseMatrix")b = 
> matrix(c(-30,40,-10),nrow=3,ncol=1)
> #solve for x, Error in LU.dgC(a) : cs_lu(A) failed: near-singular A (or out 
> of memory)solve(A,b,sparse=TRUE,tol=.Machine$double.eps)
> #one x that happens to solve Ax = bx = matrix(c(-10,10,0),nrow=3,ncol=1)A %*% 
> x
> #Error in lsfit(A, b) : only 3 cases, but 4 variableslsfit(A,b)#solves the 
> system, but fails belowsolve(qr(A, LAPACK=TRUE),b)#Error in qr.solve(A, b) : 
> singular matrix 'a' in solveqr.solve(A,b)
> #matrices used in my actual problem (see attached files)A = readMM("A.txt")b 
> = readMM("b.txt")
> #Error in as(x, "matrix")[i, , drop = drop] : subscript out of 
> boundssolve(qr(A, LAPACK=TRUE),b)

Your code is a mess. 

A singular square system of linear equations has an infinity of solutions if a 
solution exists at all.
How that works you can find here: 
https://en.wikipedia.org/wiki/System_of_linear_equations
in the section "Matrix solutions".

For your simple example you can do it like this:

library(MASS)
Ag <- ginv(A)   # pseudoinverse

xb <- Ag %*% b # minimum norm solution

Aw <- diag(nrow=nrow(Ag)) - Ag %*% A  # see the Wikipedia page
w <- runif(3)
z <- xb + Aw %*% w
A %*% z - b

N <- Null(t(A))  # null space of A;  see the help for Null in package MASS
A %*% N
A %*% (xb + 2 * N) - b

For sparse systems you will have to approach this differently; I have no 
experience with that.

Berend

______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

Reply via email to