I have a sparse matrix of type dgTMatrix
built with the Matrix package in R and I'm trying to solve for the vector x
in Ax=b
using QR decomposition, but it is not working correctly.
For example below I am solving for a random A, b
and you can see the method works with the norm of Ax-b = 0
A <- matrix(rnorm(9),ncol=3)
decomp <- qr(A)
b <- rnorm(3)
x <- qr.coef(decomp,b)
(A %*% x - b) %>% norm
[1] 3.885781e-16
My A
and b
are like this:
A: 173700 x 173700 sparse Matrix of class "dgTMatrix"```
b: 173700 x 1 sparse Matrix of class "dgCMatrix"```
When I take a QR decomposition of A
I get the following:
qr_decomp <- qr(A, LAPACK = TRUE, tol = 1e-10)
qr_decomp
'MatrixFactorization' of Formal class 'sparseQR' [package "Matrix"] with 6 slots
..@ V :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. ..@ i : int [1:13350377] 0 173283 173284 173285 173286 1 2 3 2 3 ...
.. .. ..@ p : int [1:173701] 0 5 8 11 19 22 25 30 33 36 ...
.. .. ..@ Dim : int [1:2] 173700 173700
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : NULL
.. .. ..@ x : num [1:13350377] -9.27e 01 2.35e-03 -1.28e 02 2.35e-03 6.40e 01 ...
.. .. ..@ factors : list()
..@ beta: num [1:173700] 6.88e-05 7.80e-14 1.24e-06 3.88e-05 7.80e-14 ...
..@ p : int [1:173700] 34977 38066 38838 39610 38067 38839 34980 38069 38841 39613 ...
..@ R :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. ..@ i : int [1:24075359] 0 1 1 2 0 1 2 3 4 4 ...
.. .. ..@ p : int [1:173701] 0 1 2 4 8 9 11 12 13 15 ...
.. .. ..@ Dim : int [1:2] 173700 173700
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : NULL
.. .. .. ..$ : NULL
.. .. ..@ x : num [1:24075359] 1.57e 02 3.58e 06 -1.26e 03 3.58e 06 1.92e-03 ...
.. .. ..@ factors : list()
..@ q : int [1:173700] 35749 38838 38066 36522 38839 38067 35752 38841 38069 36525 ...
..@ Dim : int [1:2] 173700 173700
According to the docs for sparseQR
class,
For a sparse m×n (“long”: m≥n) rectangular matrix A, the sparse QR decomposition is either of the form PA=QR with a (row) permutation matrix P, (encoded in the p slot of the result) if the q slot is of length 0, or of the form PAP∗=QR with an extra (column) permutation matrix P∗ (encoded in the q slot).
Since qr_decomp@q
exists, per the docs, this must be of the form PAP∗=QR
Anyway, the same procedure as above is not working in this case:
decomp <- qr(A)
x <- qr.coef(decomp, b)
(A %*% x - b) %>% norm
[1] 3.540814e 24
The above value is really far from zero.
Inspecting matrix A everything looks okay to me:
!(A@x %>% is.finite) %>% any
[1] FALSE
!(b@x %>% is.finite) %>% any
[1] FALSE
> A[1:10, 1:10]
10 x 10 sparse Matrix of class "dgTMatrix"
[1,] -5.366088e 08 5.301600e-02 . . . . . . . .
[2,] 4.418000e-02 -5.366088e 08 4.418000e-02 . . . . . . .
[3,] . 4.123467e-02 -5.366088e 08 4.123467e-02 . . . . . .
[4,] . . 3.976200e-02 -5.366088e 08 3.976200e-02 . . . . .
[5,] . . . 3.887840e-02 -5.366088e 08 3.887840e-02 . . . .
[6,] . . . . 3.828933e-02 -5.366088e 08 3.828933e-02 . . .
[7,] . . . . . 3.786857e-02 -5.366088e 08 3.786857e-02 . .
[8,] . . . . . . 2.650800e-02 -5.366088e 08 1.251767e-02 .
[9,] . . . . . . . 9.719600e-03 -5.366088e 08 9.719600e-03
[10,] . . . . . . . . 9.572333e-03 -5.366088e 08
However, there are large values present in the matrices:
max(b)
[1] 3.978441e 22
max(A)
[1] 3.979517e 22
min(b)
[1] 0
min(A)
[1] -7.958754e 22
Could it be the large values messing things up with some kind of round off error? Or am I missing some key info about QR not converging in some case?
CodePudding user response:
tl;dr it does appear that large values will lead to problems; scaling the inputs appears to help.
An example that isn't problematic, even though I used Cauchy-distributed samples to fill in the matrix (which leads to a wide range of values):
set.seed(101)
library(Matrix)
## d <- 173700
d <- 1e3
## n <- 1e6
n <- 1e5
randfun <- rcauchy
As <- sparseMatrix(i=integer(0),
j=integer(0),
dims = c(d, d), repr = "T")
As[cbind(sample(d, size = n, replace = TRUE),
sample(d, size = n, replace = TRUE))] <- randfun(n)
bs <- sparseMatrix(i = 1:d, j = rep(1, d), x = randfun(d))
qr_decomp <- qr(As, LAPACK = TRUE, tol = 1e-10)
decomp <- qr(As)
xs <- qr.coef(decomp, bs)
(As %*% xs - bs) |> norm() ## 1.19e-09
If I use larger values for d
(matrix/vector dimensions) and n
(number of non-zero elements), either the computation takes a long time or I get a segmentation fault.
However, if I set randfun <- function(n) rcauchy(n, scale = 1e20)
in the example above I do get large deviations.
Scaling the inputs seems to solve the problem, although I don't know if it will work for what you want to do ...
As2 <- As/1e20
bs2 <- bs/1e20
decomp2 <- qr(As2)
xs2 <- qr.coef(decomp2, bs2)
(As2 %*% xs2 - bs2) |> norm()