Trung tâm cá nhân
sách gpt4 ăn đã đi

c++ - 如何构造这个分块三对角(稀疏)矩阵？

In lại 作者：太空狗 更新时间：2023-10-29 21:33:10 31 4
mua khóa gpt4 giày nike

在 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

TRONG MỘT、B、C 是 5x5 矩阵，0 是 5x5 零矩阵。

实际上，我使用的矩阵是成百上千行和列。在 R 中，我知道可以使用 rbind Và cbind 但这是一个有点乏味且昂贵的解决方案。

更新:如何使用这个矩阵

令上述矩阵为H。给定两个 vector x Và S，我需要计算 H %*% x + s = y。

câu trả lời hay nhất

你如何使用这个矩阵其实更重要。在许多情况下，后续计算不需要显式矩阵构造。此问答可能与您无关:How to build & store this large lower triangular matrix for matrix-vector multiplication? , 但完美地说明了我的观点。

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

矩阵仅用于矩阵 vector 乘法？我们绝对可以跳过形成这个矩阵，因为乘法只是 rbind(B, A, C) Và x 之间的滚动矩阵 vector 乘法。

## `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)

为了验证上面的Và是否被正确计算，我们需要显式构造H。构建方式有很多种。

方法一:使用分块对角(稀疏)矩阵

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 是正确的。

方法二:直接填写

此方法基于以下观察:

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++ .

关于c++ - 如何构造这个分块三对角(稀疏)矩阵？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/52452554/

31 4 0
Đề xuất bài viết: python - 如何处理具有多个窗口的林间空地项目
Đề xuất bài viết: c# - 为什么 KeyPress 的消息框显示在 KeyDown 之前？
Đề xuất bài viết: c++ - Python C-API 对象分配
Đề xuất bài viết: c# - 启动 Windows 10 商店应用程序
太空狗
Hồ sơ

Tôi là một lập trình viên xuất sắc, rất giỏi!

Nhận phiếu giảm giá taxi Didi miễn phí
Phiếu giảm giá taxi Didi
Các bài viết phổ biến trên toàn bộ trang web
Chứng chỉ ICP Bắc Kinh số 000000
Hợp tác quảng cáo: 1813099741@qq.com 6ren.com
Xem sitemap của VNExpress