Как изменить строки матрицы в зависимости от строки (Создание матрицы Хаара) - PullRequest
Новая
       8

Как изменить строки матрицы в зависимости от строки (Создание матрицы Хаара)

0 голосов
/ 01 декабря 2018

Я пишу пакет, который требует создания очень больших матриц Хаара ( $ 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 {bmatrix} $$

(обратите внимание, что размер всегда является степенью 2)

Я написал функцию, которая будет строить матрицу Хаара произвольного размера (2 ^ J) с использованием продуктов Кронекера.

$$ 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 ((строка - 1) - 2 ^ {\ lfloor \ log_2(строка - 1) \ rfloor} \ right) $$

Таким образом,

  • positiveOneStart = width * (row - 1 - 2 ^ этаж (журнал (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

Может кто-нибудь помочь мне понять, что не так?Я проверил эту штуку до смерти, и я не могу понять, где ошибка.

...