How to construct this block tridiagonal (sparse) matrix?

Question

Is there a quick way to fill a (sparse) matrix as such in either R or C++:

A, B, 0, 0, 0
C, A, B, 0, 0
0, C, A, B, 0
0, 0, C, A, B
0, 0, 0, C, A

Where A, B, C are 5x5 matrices and 0 is a 5x5 matrix of zeros.

In reality, the matrices I use are hundreds to thousands of rows by columns. In R, I know that rbind and cbind could be used but this is a somewhat tedious and expensive solution.

update: How this matrix is to be used

Let the above matrix be H. Given two vectors x and s, I need to compute H %*% x + s = y.

depends... how do you save your matrices, which container, sparse or not. A bit more information would be useful... (c++) — , Sep 22 '18 at 01:01
In c++ I normally use Eigen. Most of the time the matrix is sparse, but sometimes dense when added to another dense matrix. — Tiffany, Sep 22 '18 at 01:54
If we call the above matrix M, have a vector x, and a vector s: M*x + s = y, which is computed each timestep where y becomes x in the next iteration. — Tiffany, Sep 22 '18 at 02:44
I normally code in R, and I'm building a simple small grid of this model for testing in R. Then I normally convert it into c++ for speed, but I'm willing to skip the R part if absolutely necessary. I'm not nearly as proficient in c++, so debugging and working through the math is easier for me in R. — Tiffany, Sep 22 '18 at 02:45
Correct. There's an implicit version of the model which would require inverting or solving the system (solving for x and knowing y), but for this version only a matrix vector multiplication plus addition is required. — Tiffany, Sep 22 '18 at 03:24

Zheyuan Li · Accepted Answer · 2018-09-22T04:19:13.587

How you use this matrix is actually more important. In many cases, explicit matrix construction is not necessary for subsequent computations. This Q & A may not be relevant to you: How to build & store this large lower triangular matrix for matrix-vector multiplication?, but is perfect to illustrate my point.

Let the above matrix be H. Given two vectors x and s, I need to compute H %*% x + s = y.

The matrix is only used in a matrix-vector multiplication? We can definitely skip forming this matrix, as the multiplication is just a rolling matrix-vector multiplication between rbind(B, A, C) and x.

## `nA` is the number of `A`-blocks on the main diagonal of `H`
MatVecMul <- function (A, B, C, nA, x, s) {
  ## input validation
  if (diff(dim(A))) stop("A is not a square matrix")
  if (diff(dim(B))) stop("B is not a square matrix")
  if (diff(dim(C))) stop("C is not a square matrix")
  if (dim(A)[1] != dim(B)[1]) stop("A and B does not have the same dimension")
  if (dim(A)[1] != dim(C)[1]) stop("A and C does not have the same dimension")
  if (length(x) != nA * M) stop("dimension dismatch between matrix and vector")
  if (length(x) %% length(s)) stop("length of 'x' does not divide length of 's'")
  ## initialization
  y <- numeric(length(x))
  ##########################
  # compute `y <- H %*% x` #
  ##########################
  ## first block column contains `rbind(A, C)`
  M <- dim(A)[1]
  ind_x <- 1:M
  y[1:(2 * M)] <- rbind(A, C) %*% x[ind_x]
  ind_x <- ind_x + M
  ## middle (nA - 2) block columns contain `rbind(B, A, C)`
  BAC <- rbind(B, A, C)
  ind_y <- 1:(3 * M)
  i <- 0
  while (i < (nA - 2)) {
    y[ind_y] <- y[ind_y] + BAC %*% x[ind_x]
    ind_x <- ind_x + M
    ind_y <- ind_y + M
    i <- i + 1
    }
  ## final block column contains `rbind(A, C)`
  ind_y <- ind_y[1:(2 * M)]
  y[ind_y] <- y[ind_y] + rbind(B, A) %*% x[ind_x]
  ## compute `y + s` and return
  y + s
  }

Here is a reproducible example.

set.seed(0)
M <- 5  ## dim of basic block
A <- matrix(runif(M * M), M)
B <- matrix(runif(M * M), M)
C <- matrix(runif(M * M), M)
nA <- 5
x <- runif(25)
s <- runif(25)

y <- MatVecMul(A, B, C, nA, x, s)

To verify that the above y is correctly computed, we need to explicitly construct H. There are many ways for construction.

method 1: use block diagonal (sparse) matrix

N <- nA * M  ## dimension of the final square matrix

library(Matrix)

## construct 3 block diagonal matrices
H1 <- bdiag(rep.int(list(A), nA))
H2 <- bdiag(rep.int(list(B), nA - 1))
H3 <- bdiag(rep.int(list(C), nA - 1))

## augment H2 and H3, then add them together with H1
H <- H1 +
     rbind(cbind(Matrix(0, nrow(H2), M), H2), Matrix(0, M, N)) + 
     cbind(rbind(Matrix(0, M, ncol(H3)), H3), Matrix(0, N, M))

## verification
range((H %*% x)@x + s - y)
#[1] -8.881784e-16  8.881784e-16

We see that MatVecMul is correct.

method 2: direct fill-in

This method is based on the following observation:

B
-------------
A  B
C  A  B
   C  A  B
      C  A  B
         C  A
-------------
            C

It is easy to construct the rectangular matrix first, then subset the square matrix in the middle.

BAC <- rbind(B, A, C)

nA <- 5  ## number of basic block
N <- nA * M  ## dimension of the final square matrix
NR <- N + 2 * M  ## leading dimension of the rectangular matrix

## 1D index for the leading B-A-C block
BAC_ind1D <- c(outer(1:nrow(BAC), seq(from = 0, by = NR, length = M), "+"))
## 1D index for none-zero elements in the rectangular matrix
fill_ind1D <- outer(BAC_ind1D, seq(from = 0, by = M * (NR + 1), length = nA), "+")
## 2D index for none-zero elements in the rectangular matrix
fill_ind2D <- arrayInd(fill_ind1D, c(NR, N))

## construct "dgCMatrix" sparse matrix
library(Matrix)
Hsparse <- sparseMatrix(i = fill_ind2D[, 1], j = fill_ind2D[, 2], x = BAC)
Hsparse <- Hsparse[(M+1):(N+M), ]

## construct dense matrix
Hdense <- matrix(0, NR, N)
Hdense[fill_ind2D] <- BAC
Hdense <- Hdense[(M+1):(N+M), ]

## verification
range((Hsparse %*% x)@x + s - y)
#[1] -8.881784e-16  8.881784e-16

range(base::c(Hdense %*% x) + s - y)
#[1] -8.881784e-16  8.881784e-16

Once again, we see that MatVecMul is correct.

implementing MatVecMul with Rcpp

It is very easy to transform the R function MatVecMul to an Rcpp function. I would leave this task to you as you have used c++.

How to construct this block tridiagonal (sparse) matrix?

1 Answers1