[R] Solving sparse, singular systems of equations

[Jeff Newmiller]

The usual culprit in messy code is posting in HTML format. That usually leads to stripping of the formatting by the mailing list and a notice that that occurred, but I don't see that warning here. I still think posting plain text format would fix the problem. 
-- 
Sent from my phone. Please excuse my brevity.

On April 20, 2016 11:51:40 AM PDT, A A via R-help <r-help at r-project.org> wrote:
>Thanks for the help. Sorry, I am not sure why it looks like that in the
>mailing list - it looks much more neat on my end (see attached file). 
>
>On Wednesday, April 20, 2016 2:01 PM, Berend Hasselman <bhh at xs4all.nl>
>wrote:
> 
>
> 
>> On 20 Apr 2016, at 13:22, A A via R-help <r-help at 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 at 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.

	[[alternative HTML version deleted]]