I have a sparse quadratic diagonal matrix H of dimension (K*M) x (K*M). The goal is to replace the lower off-diagonal blocks of the matrix with a known quadratic matrix A with size K x K:
# dimensions of the final matrix are (K*M) x (K*M)
K = 5
M = 3
# initiate the matrix as diagonal
H = Matrix::Diagonal(K*M)
# convert because off-diagonal replacement does not work for sparse diagonal class
H = as(H, "dgCMatrix")
# matrix that will replace the lower off-diagonal blocks of size K x K
A = matrix(runif(K^2), K, K)
The following method uses an auxiliary vector that indexes the off-diagonal blocks of appropriate size and loops over each block:
# construct vector that is used to index the lower off-diagonal blocks of size K x K
indices = seq(1, K * (1 + M), K)
# replace all M - 1 lower off-diagonal blocks
for(ii in 1:(M-1)){
H[indices[ii + 1]:(indices[ii + 2]-1), indices[ii]:(indices[ii + 1]-1)] = A
}
Unfortunately, indexing and replacing these blocks becomes extremely slow when the dimensionality of the matrix grows. Are there faster methods available that would work reasonably well in settings with, say, K=8000 and M=20 ? It does not necessarily need to be within the Matrix package setting as long as the result can be back-converted to a Matrix class.
The Kronecker product of a sparse banded matrix and the dense block should be "fast", but you'll want to do your own benchmarks ...
library(Matrix)
set.seed(0)
K <- 5L
M <- 3L
(X <- bandSparse(M, M, -1L, list(rep.int(1, M - 1L))))
## 3 x 3 sparse Matrix of class "dtCMatrix"
##
## [1,] . . .
## [2,] 1 . .
## [3,] . 1 .
(Y <- matrix(round(runif(K * K), 2L), K, K))
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.90 0.20 0.06 0.77 0.78
## [2,] 0.27 0.90 0.21 0.50 0.93
## [3,] 0.37 0.94 0.18 0.72 0.21
## [4,] 0.57 0.66 0.69 0.99 0.65
## [5,] 0.91 0.63 0.38 0.38 0.13
XY <- kronecker(X, Y)
diag(XY) <- 1
XY
## 15 x 15 sparse Matrix of class "dgCMatrix"
##
## [1,] 1.00 . . . . . . . . . . . . . .
## [2,] . 1.00 . . . . . . . . . . . . .
## [3,] . . 1.00 . . . . . . . . . . . .
## [4,] . . . 1.00 . . . . . . . . . . .
## [5,] . . . . 1.00 . . . . . . . . . .
## [6,] 0.90 0.20 0.06 0.77 0.78 1.00 . . . . . . . . .
## [7,] 0.27 0.90 0.21 0.50 0.93 . 1.00 . . . . . . . .
## [8,] 0.37 0.94 0.18 0.72 0.21 . . 1.00 . . . . . . .
## [9,] 0.57 0.66 0.69 0.99 0.65 . . . 1.00 . . . . . .
## [10,] 0.91 0.63 0.38 0.38 0.13 . . . . 1.00 . . . . .
## [11,] . . . . . 0.90 0.20 0.06 0.77 0.78 1 . . . .
## [12,] . . . . . 0.27 0.90 0.21 0.50 0.93 . 1 . . .
## [13,] . . . . . 0.37 0.94 0.18 0.72 0.21 . . 1 . .
## [14,] . . . . . 0.57 0.66 0.69 0.99 0.65 . . . 1 .
## [15,] . . . . . 0.91 0.63 0.38 0.38 0.13 . . . . 1