2

I'm trying to decompose a square matrix with the R package irlba but am getting the following message:

"Error in V[, 1:(k + dim(F)[2])] <- cbind(V[, 1:(dim(Bsvd$v)[1]), drop = FALSE] %*% : number of items to replace is not a multiple of replacement length"

Decomposing the matrix with svd or eigenwork, however. I found this post elsewhere alluding to the same problem (without a response).

I would be grateful to anyone that can help me figure this out!

Code:

library(irlba)

C.i <- structure(c(0.107821513562202, 0.0629415996952743, -0.0346069282137902, 
-0.0410208578535759, 0.0629415996952743, 0.136205395050303, -0.00378166611862939, 
0.00237156895349009, -0.0346069282137902, -0.00378166611862939, 
0.0122114763151942, -0.00635448918784993, -0.0410208578535759, 
0.00237156895349009, -0.00635448918784993, 0.0431193044924), .Dim = c(4L, 4L), 
.Dimnames = list(c("Sepal.Length", "Sepal.Width", "Petal.Length",
"Petal.Width"), c("Sepal.Length", "Sepal.Width", "Petal.Length", 
"Petal.Width")))

irlba(C.i, nu=1, nv=1)

# These work
svd(C.i)
eigen(C.i)
Community
  • 1
  • 1
Marc in the box
  • 11,769
  • 4
  • 47
  • 97
  • I usually dig the source and try to throw out everything that's unrelated until I have a minimum breaking code. It is usually enough to diagnose if that's a bug or e.g. incorrect params. See [this](http://stackoverflow.com/questions/25476726/why-doesnt-solve-qp-and-portfolio-optim-generate-identical-results) question, may give a methodological hint. If you manage to narrow it down at least partially, do not hesitate to post, I might pick it up from there when I have some time. – tonytonov Sep 01 '14 at 09:35
  • I thought this might be an edge case having to do with the fact that `adjust+max(nu,nv)` is equal to the dimension of the matrix, but that doesn't seem to be true: `irlba(C.i, nu=1,nv=1,adjust=2)` doesn't work either. However, `irlba(matrix(rnorm(25),5),nv=1,nu=1)` works and `irlba(matrix(rnorm(16),4),nv=1,nu=1)` fails -- in case that helps ... – Ben Bolker Sep 07 '14 at 22:34
  • @BenBolker - Yes - for some reason, there are only some matrices that result in these types of problems. I think it may have do do with the situations where the matrix is not positive definite, but I'm not sure about that. – Marc in the box Sep 08 '14 at 09:11
  • @tonytonov - I had limited luck going through the original function line by line to try and find the point of error. I have since contacted the author of the package and he is also taking a look at this error. Thanks for the advise though. – Marc in the box Sep 08 '14 at 09:12

2 Answers2

1

With irlba(C.i, nu=0, nv=0) I get:

$d
[1] 0.1938809

$u
           [,1]
[1,]  0.6595482
[2,]  0.7216141
[3,] -0.1355695
[4,] -0.1609040

$v
           [,1]
[1,]  0.6590202
[2,]  0.7214017
[3,] -0.1432499
[4,] -0.1573256

$iter
[1] 1

$mprod
[1] 6

Is this what you were looking for?

(Sorry for asking a question in an answer, but I don't have enough points to comment on the question. :| )

Yay295
  • 1,628
  • 3
  • 17
  • 29
  • This does indeed work, but I don't know why. I've tried it out on one example and it gives the same result as `nu=1`and `nv=1`: `set.seed(1); x <- matrix(rnorm(25),5); irlba(x,nv=1,nu=1); irlba(x,nv=0,nu=0)`. Do you know what the function is doing under these settings=? – Marc in the box Sep 08 '14 at 10:53
  • I had a look at the code, but I'm not familiar enough with the functions to really know what's going on. I think it may just be a rare bug. As to why it works with 0, I think that may be another bug. If you'll notice in the test you did there ^, u and v are the same, but iter and mprod are not. While this does seem to work here, it may not work for another matrix. In fact, I just found a random matrix whose results are all of the wrong sign, so be careful. – Yay295 Sep 09 '14 at 07:28
  • Yes, I believe you are right that this not a good overall fix. I tried the `nu=0; nv=0` trick with a much larger matrix and the function continued to calculate for a much longer time than it should have if only a single vector was to be computed. I eventually had to break the calculation. I will keep the question open and hope for a more canonical answer. Thanks for your help. – Marc in the box Sep 09 '14 at 12:40
0

I have received a response from the package author regarding the error - It was indeed a bug. This problem has been remedied with an updated version of the irlba package. The updated package will eventually be submitted to CRAN, but for the meantime you can install via GitHub:

remove.packages("irlba")
library("devtools")
install_github("IRL","bwlewis",quick=TRUE)

library("irlba")

C.i <- structure(c(0.107821513562202, 0.0629415996952743, -0.0346069282137902, 
-0.0410208578535759, 0.0629415996952743, 0.136205395050303, -0.00378166611862939, 
0.00237156895349009, -0.0346069282137902, -0.00378166611862939, 
0.0122114763151942, -0.00635448918784993, -0.0410208578535759, 
0.00237156895349009, -0.00635448918784993, 0.0431193044924), .Dim = c(4L, 4L), 
.Dimnames = list(c("Sepal.Length", "Sepal.Width", "Petal.Length",
"Petal.Width"), c("Sepal.Length", "Sepal.Width", "Petal.Length", 
"Petal.Width")))

irlba(C.i, nu=1, nv=1)

#$d
#[1] 0.1938898
#
#$u
#           [,1]
#[1,]  0.6593389
#[2,]  0.7216000
#[3,] -0.1349355
#[4,] -0.1623519
#
#$v
#           [,1]
#[1,]  0.6593384
#[2,]  0.7216001
#[3,] -0.1349364
#[4,] -0.1623526
#
#$iter
#[1] 2
#
#$mprod
#[1] 10
Marc in the box
  • 11,769
  • 4
  • 47
  • 97