如何为矩阵向量乘法构建和存储这个大的下三角矩阵？

Question

我需要创建一个具有特殊结构的下三角矩阵，然后进行矩阵向量乘法。

矩阵由值 k 参数化。它的主对角线是k ^ 0的向量，即1；第一个次对角线是 k ^ 1 的向量，第 i 个次对角线包含 k ^ i.

这是一个 5 x 5 的示例 k = 0.9:

structure(c(1, 0.9, 0.81, 0.729, 0.6561, 0, 1, 0.9, 0.81, 0.729, 
0, 0, 1, 0.9, 0.81, 0, 0, 0, 1, 0.9, 0, 0, 0, 0, 1), .Dim = c(5L, 5L))
#       [,1]  [,2] [,3] [,4] [,5]
#[1,] 1.0000 0.000 0.00  0.0    0
#[2,] 0.9000 1.000 0.00  0.0    0
#[3,] 0.8100 0.900 1.00  0.0    0
#[4,] 0.7290 0.810 0.90  1.0    0
#[5,] 0.6561 0.729 0.81  0.9    1

我需要构造一个100,000 x 100,000这样大的矩阵并用它来计算。为此，我需要最有效的存储方法。有任何想法吗？

Answer 1

试试这个：

k <- 0.9
n <- 5

d <- diag(n)
replace(k ^ (row(d) - col(d)), upper.tri(d), 0)

给予：

       [,1]  [,2] [,3] [,4] [,5]
[1,] 1.0000 0.000 0.00  0.0    0
[2,] 0.9000 1.000 0.00  0.0    0
[3,] 0.8100 0.900 1.00  0.0    0
[4,] 0.7290 0.810 0.90  1.0    0
[5,] 0.6561 0.729 0.81  0.9    1

Answer 2

您并不总是需要明确地形成矩阵来进行矩阵-向量或矩阵-矩阵乘法。例如，没有人真正形成对角矩阵并将其用于此类计算。

你的矩阵和对角矩阵没有本质区别。

所以你将操作简化为一系列向量加法。这是一个简单的 R 级实现。

MatVecMul <- function (y, k) {
  n <- length(y)
  z <- numeric(n)
  for (i in 1:n) z[i:n] <- z[i:n] + k ^ (i - 1) * y[1:(n - i + 1)]
  z
  }

与直接矩阵构造和计算的比较。

d <- structure(c(1, 0.9, 0.81, 0.729, 0.6561, 0, 1, 0.9, 0.81, 0.729, 
0, 0, 1, 0.9, 0.81, 0, 0, 0, 1, 0.9, 0, 0, 0, 0, 1), .Dim = c(5L, 5L))
set.seed(0); y <- runif(5)
c(d %*% y)
#[1] 0.8966972 1.0725361 1.3374064 1.7765191 2.5070750

MatVecMul(y, 0.9)
#[1] 0.8966972 1.0725361 1.3374064 1.7765191 2.5070750

可以用 Rcpp 轻松替换 R 级 "for" 循环。

library(Rcpp)
cppFunction("NumericVector MatVecMul_cpp (NumericVector y, double k) {
  int n = y.size();
  NumericVector z(n);
  int i; double *p1, *p2, *end = &z[n];
  double tmp = 1.0;
  for (i = 0; i < n; i++) {
    for (p1 = &z[i], p2 = &y[0]; p1 < end; p1++, p2++) *p1 += tmp * (*p2);
    tmp *= k;
    }
  return z;
  }")

MatVecMul_cpp(y, 0.9)
#[1] 0.8966972 1.0725361 1.3374064 1.7765191 2.5070750

让我们有一个基准。

v <- runif(1e4)
system.time(MatVecMul(y, 0.9))
#   user  system elapsed 
#  3.196   0.000   3.198 
system.time(MatVecMul_cpp(y, 0.9))
#   user  system elapsed 
#  0.840   0.000   0.841

但要注意一点：注意机器精度。一旦 k ^ (i - 1) 变得太小，您可能会在加法过程中丢失所有有效数字。参见 . In this example with k = 0.9, there is k ^ 400 = 5e-19. So even though the full matrix is 10000 x 10000, it is numerically banded 比下三角。这意味着我们实际上可以更早地终止循环。但是我不会实现这个。

如何为矩阵向量乘法构建和存储这个大的下三角矩阵？

How to build & store this large lower triangular matrix for matrix-vector multiplication?

performance

r

matrix

out-of-memory

matrix-multiplication