如何构建这个块三对角（稀疏）矩阵？

Question

在R或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

其中A、B和C是5x5矩阵，而0是一个5x5的零矩阵。

实际上，我使用的矩阵可能有数百到数千行和列。在R中，我知道可以使用rbind和cbind，但这是一种有点繁琐和昂贵的解决方案。

---

更新：如何使用此矩阵

设上述矩阵为H。给定两个向量x和s，我需要计算H %*% x + s = y。

- Tiffany

---

取决于你如何保存矩阵，使用的容器是稠密还是稀疏。更多信息会更有用...(c++) - user8705939

在C++中，我通常使用Eigen。大多数情况下，矩阵是稀疏的，但有时会与另一个密集矩阵相加而变得密集。 - Tiffany

如果我们将上述矩阵称为M，有一个向量x和一个向量s：M * x + s = y，在每个时间步骤计算y成为下一次迭代中的x。 - Tiffany

我通常使用R进行编码，并正在R中构建一个简单的小模型网格以进行测试。然后，我通常将其转换为C++以提高速度，但如果绝对必要，我愿意跳过R部分。我在C++方面的熟练程度远不如R，因此在R中进行调试和数学计算更容易。 - Tiffany

正确。模型有一个隐式版本，需要求逆或解决系统（解出x并知道y），但对于这个版本，只需要进行矩阵向量乘法加上加法即可。 - Tiffany

Answer 1

这个矩阵的使用方式实际上更为重要。在许多情况下，后续计算并不需要显式地构建矩阵。这个问答可能与您无关：如何构建和存储大型下三角矩阵以进行矩阵向量乘法？，但它很适合说明我的观点。

引用： 给定两个向量x和s，我需要计算H % *% x + s = y。

这个矩阵只用于矩阵向量乘法吗？我们肯定可以跳过形成这个矩阵，因为乘法只是rbind（B，A，C）和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
  }

这是一个可重现的例子。

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)

---

为了验证上述的y计算正确，我们需要明确构造H。有很多构造方法。

方法1：使用块对角线（稀疏）矩阵。

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

我们看到MatVecMul是正确的。

方法2：直接填充

该方法基于以下观察结果：

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

很容易先构建矩形矩阵，然后从中间子集出正方形矩阵。

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

再次看到，我们发现MatVecMul是正确的。

---

使用Rcpp实现MatVecMul

将R函数MatVecMul转换为Rcpp函数非常容易。我会把这个任务交给你，因为你已经使用过C++了。

- Zheyuan Li