如何建造和建造存储此较大的下三角矩阵以进行矩阵矢量乘法? [英] How to build & store this large lower triangular matrix for matrix-vector multiplication?

问题描述

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

I need to create a lower triangular matrix with a special structure then do a matrix-vector multiplication.

矩阵由值k设置参数.它的主对角线是k ^ 0的矢量，即1；第一个对角线是k ^ 1的向量，而第i个子对角线保持k ^ i.

The matrix is parameterized by a value k. It main diagonal is a vector of k ^ 0, i.e. 1; the first sub-diagonal is a vector of k ^ 1, and the i-th sub-diagonal holds k ^ i.

这是一个的5 x 5示例:

Here is a 5 x 5 example with 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一样大的矩阵，并将其用于计算.我需要最有效的存储方法.有任何想法吗?

I need to construct such a matrix as large as 100,000 x 100,000 and use it for computation. I need the most efficient storage method for this. Any ideas?

推荐答案

您不必总是显式地形成矩阵来进行矩阵矢量或矩阵矩阵乘法.例如，没有人真正形成对角矩阵并将其用于这种计算.

You don't always need to explicitly form a matrix to do a matrix-vector or matrix-matrix multiplication. For example, no one really forms a diagonal matrix and use it for such computations.

矩阵和对角矩阵之间没有实质性差异.

There is no substantial difference between your matrix and a diagonal matrix.

因此，您可以将操作简化为一系列矢量加法.这是一个简单的R级实现.

So you reduce the operation to a series of vector addition. Here is a trivial R-level implementation.

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
  }

与直接矩阵构建和计算的比较.

A comparison with direct matrix construction and computation.

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"循环.

Can replace the R-level "for" loop easily with Rcpp.

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)变得太小，在加法过程中您可能会丢失所有有效数字.参见 R:使用`(1 + 1/n)^ n`近似`e = exp(1)`会导致荒谬的结果，当`n `很大.在带有k = 0.9的示例中，存在k ^ 400 = 5e-19.因此，即使完整的矩阵为10000 x 10000，其数字也比下三角的带状.这意味着我们实际上可以更早地终止循环.但是我不会实现这一点.

One caution though: be aware of machine precision. As soon as k ^ (i - 1) becomes too small, you may lose all significant digits during addition. See R: approximating `e = exp(1)` using `(1 + 1 / n) ^ n` gives absurd result when `n` is large. 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 than lower triangular. This means that we can actually terminate the loop earlier. But I will not implement this.

这篇关于如何建造和建造存储此较大的下三角矩阵以进行矩阵矢量乘法?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！