r - 如何根据行修改矩阵的行（制作 Haar 矩阵）

Question

我正在编写一个需要制作非常大的 Haar 矩阵的包（$2^{28} \约 270$百万行和列）。例如：

$$ H_2 = \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} $$

$$ H_4 = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & 0 & 0\\ 0 & 0 & 1 & -1 \end{矩阵} $$

（请注意，维度始终是 2 的幂）

我编写了一个函数，该函数将使用 Kroneker 产品构建任意大小 (2^J) 的 Haar 矩阵。

$$ H_{2N} = \begin{bmatrix} H_{N} \otimes [1, 1] \\ I_{N} \otimes [1, -1] \end{bmatrix} $$

然而，这种算法的迭代性质在内存方面非常浪费，并且在R.

THANKFULLY，特定行的哪些值应该是正数或负数的表达式可以以封闭形式给出，如下所示：

---

• 对于给定值$J$，初始化一个维度为零的矩阵$2^J \times 2^J$

• 将第一行设置为所有。

• 对于所有其他行，width任何给定行中的非零数字序列的 是
  $$2^{J - \lfloor \log_2(row - 1) \rfloor} $$

• 上述组的start位置是
  $$\text{width} * \left((row - 1) - 2^{\lfloor \log_2(row - 1) \rfloor}\right)$$

因此，

• positiveOneStart= width* ( row- 1 - 2^floor(log( row- 1, 2))) + 1
• finalOnePosition= ( start+ width/2 - 1)
• negativeOneStart= finalOnePosition+ 1
• finalNegativeOnePosition= positiveOneStart+ width- 1

---

我写了以下两个函数：

rowFill <- function(row, H){

  if(row == 1){
    H[row, ] <<- 1
  }else{
    (log2row <- 2^floor(log(row - 1, 2)))
    (width <- dimension / log2row %>%
      round())
    (width <- suppressWarnings(as.integer(width)))

    (start = width * (row - 1 - 2^floor(log(row - 1, 2))) + 1 %>%
      round())
    (startOnes <- as.integer(start))

    (endOnes <- startOnes + width/2 -1)

    (startNegOnes <-  endOnes + 1)

    (endNegOnes <- startOnes + width - 1)

    H[row, startOne:endOnes] <<- (1)
    H[row, startNegOnes:endNegOnes] <<- (-1)
  }
}

当我在循环中使用上面的想法来填充行时，它似乎工作正常。

library(purrr)

J = 3
dimension <- 2^J 
row <- 1:(dimension)
H = Matrix::Matrix(0, nrow = dimension, ncol = dimension) %>% methods::as("dgCMatrix")

H[1,] <- 1

for(row in 2:2^J){
width <- dimension / 2^floor(log(row - 1, 2)) %>%
  round()
width <- suppressWarnings(as.integer(width))

start = width * (row - 1 - 2^floor(log(row - 1, 2))) + 1 %>%
  round()
start <- as.integer(start)

endOnes <- start + width/2 -1

startNegOnes <-  endOnes + 1

endNegOnes <- start + width - 1

H[row,start:endOnes] <- 1
H[row,startNegOnes:endNegOnes] <- -1

}

H

有用！

8 x 8 sparse Matrix of class "dgCMatrix"

[1,] 1  1  1  1  1  1  1  1
[2,] 1  1  1  1 -1 -1 -1 -1
[3,] 1  1 -1 -1  .  .  .  .
[4,] .  .  .  .  1  1 -1 -1
[5,] 1 -1  .  .  .  .  .  .
[6,] .  .  1 -1  .  .  .  .
[7,] .  .  .  .  1 -1  .  .
[8,] .  .  .  .  .  .  1 -1

但是，我想避免循环，因为我正在处理一个$2^{28}$（非常大）的矩阵，并且循环它们可能非常慢。所以我想我会尝试使用purrr包函数walk来修改 H 矩阵作为副作用：

haarFill <- function(J){
  dimension <- 2^J 
  row <- 1:(dimension)
  H = Matrix::Matrix(0, nrow = dimension, ncol = dimension) %>% methods::as("dgCMatrix")
  walk(.x = row, .f = rowFill, H)
}

但是，这并没有给我想要的结果：

8 x 8 sparse Matrix of class "dgCMatrix"

[1,] 1  1  1  1  1  1  1  1
[2,] .  .  .  1 -1 -1 -1 -1
[3,] .  1 -1 -1  1  .  .  .
[4,] .  .  .  .  1  1 -1 -1
[5,] 1 -1  1  1  1  .  .  .
[6,] .  .  1 -1  1  .  .  .
[7,] .  .  .  .  1 -1  .  .
[8,] .  .  .  .  1  1  1 -1

谁能帮我看看有什么问题？我已经对这个东西进行了单元测试，但我似乎无法弄清楚错误在哪里。

标签： ralgorithmmatrixwavelet