Написание функции факторизации QR домовладельца в коде R

Question

Написание функции факторизации QR домовладельца в коде R

Я работаю над фрагментом кода, чтобы найти QR-факторизацию матрицы в R.

X <- structure(c(0.8147, 0.9058, 0.127, 0.9134, 0.6324, 0.0975, 0.2785, 
0.5469, 0.9575, 0.9649, 0.1576, 0.9706, 0.9572, 0.4854, 0.8003
), .Dim = c(5L, 3L))


myqr <- function(A) {
  n <- nrow(A)
  p <- ncol(A)
  Q <- diag(n)
  Inp <- diag(nrow = n, ncol = p)

  for(k in c(1:ncol(A))) {
    # extract the kth column of the matrix
    col<-A[k:n,k]
    # calculation of the norm of the column in order to create the vector
    norm1<-sqrt(sum(col^2))
    # Define the sign positive if a1 > 0 (-) else a1 < 0(+)  
    sign <- ifelse(col[1] >= 0, -1, +1)  
    # Calculate of the vector a_r
    a_r <- col - sign * Inp[k:n,k] * norm1
    # beta = 2 / ||a-r||^2  
    beta <- 2 / sum(t(a_r) %*% a_r)
    # the next line of code calculates the matrix Q in every step
    Q <- Q - beta *Q %*% c(rep(0,k-1),a_r) %*% t(c(rep(0,k-1),a_r))    
    # calculates the matrix R in each step
    A[k:n,k:p] <- A[k:n,k:p] - beta * a_r %*% t(a_r) %*% A[k:n,k:p]
    }

  list(Q=Q,R=A)
  }

Но, здесь я не рассчитал в каждом шаге матрицы H что представляет отражение домохозяина, также я не рассчитал матрицу A на каждом шагу.

Как H = I - 2 v v', если я умножить на Q Я получаю

QH = Q - 2 (Qv) v'    // multiplication on the left
HQ = Q - 2 v (Q'v)'    // multiplication on the right

Теперь эти операции должны работать на каждом шагу. Однако если я рассмотрю первую матрицу H а он вторая матрица H1.... эти матрицы будут меньше, чем первая. Чтобы избежать этого, я использовал следующую строку кода:

 Q <- Q - beta * Q %*% c(rep(0,k-1),a_r) %*% t(c(rep(0,k-1),a_r))

но я не уверен, почему код работает хорошо, когда я генерирую новый вектор a_r с первым k записи нулей на каждом шагу.

3

r matrix qr-decomposition

Источник

user3483060 04 окт '16 в 10:28

1 ответ

Решение

Другие вопросы по тегам r matrix qr-decomposition

user4891738 04 окт '16 в 11:39 2016-10-04 11:39 · Accepted Answer · 2016-10-04 11:39

Я думал, что вы хотите точно такой же вывод, как возвращается qr.default, который использует компактное хранилище QR. Но потом я понял, что вы храните Q а также R факторы отдельно.

Обычно QR-факторизация только формирует R но нет Q, Далее я опишу QR-факторизацию, где формируются оба. Для тех, кому не хватает базовых знаний о факторизации QR, сначала прочтите это: lm (): что такое qraux, возвращаемое QR-декомпозицией в LINPACK / LAPACK, где в LaTeX есть аккуратные математические формулы. В дальнейшем я буду предполагать, что каждый знает, что такое отражение домохозяина и как оно вычисляется.

Процедура QR-факторизации

Во-первых, вектор рефлексии домохозяина H = I - beta * v v' (где beta вычисляется как в вашем коде), а не H = I - 2 * v v',

Затем QR-факторизация A = Q R продолжается как (Hp ... H2 H1) A = R, где Q = H1 H2 ... Hp, Вычислить Qинициализируем Q = I (единичная матрица), затем умножить Hk справа итеративно в цикле. Чтобы вычислить R, мы инициализируем R = A и умножить Hk слева итеративно в цикле.

Теперь, на k-й итерации, мы имеем обновление матрицы ранга 1 Q а также A:

Q := Q Hk = Q (I - beta v * v') = Q - (Q v) (beta v)'
A := Hk A = (I - beta v * v') A = A - (beta v) (A' v)'

v = c(rep(0, k-1), a_r), где a_r является уменьшенной ненулевой частью вектора полного отражения.

Код, который у вас есть, выполняет такое обновление с жесткой силой:

Q <- Q - beta * Q %*% c(rep(0,k-1), a_r) %*% t(c(rep(0,k-1),a_r))

Это первые колодки a_r получить полный вектор отражения и выполнить обновление ранга 1 для всей матрицы. Но на самом деле мы можем сбросить эти нули и написать (сделайте некоторую матричную алгебру, если неясно):

Q[,k:n] <- Q[,k:n] - tcrossprod(Q[, k:n] %*% a_r, beta * a_r)
A[k:n,k:p] <- A[k:n,k:p] - tcrossprod(beta * a_r, crossprod(A[k:n,k:p], a_r))

так что только часть Q а также A обновляются.

Несколько других комментариев к вашему коду

Вы использовали t() а также "%*%" много! Но почти все они могут быть заменены crossprod() или же tcrossprod(), Это исключает явное транспонирование t() и более эффективно использовать память;
Вы инициализируете другую диагональную матрицу Inp что не обязательно. Получить вектор отражения домохозяина a_rможно заменить
```
sign <- ifelse(col[1] >= 0, -1, +1)
a_r <- col - sign * Inp[k:n,k] * norm1
```
от
```
a_r <- col; a_r[1] <- a_r[1] + sign(a_r[1]) * norm1
```
где sign является базовой функцией R.

R код для QR факторизации

## QR factorization: A = Q %*% R
## if `complete = FALSE` (default), return thin `Q`, `R` factor
## if `complete = TRUE`, return full `Q`, `R` factor

myqr <- function (A, complete = FALSE) {

  n <- nrow(A)
  p <- ncol(A)
  Q <- diag(n)

  for(k in 1:p) {
    # extract the kth column of the matrix
    col <- A[k:n,k]
    # calculation of the norm of the column in order to create the vector r
    norm1 <- sqrt(drop(crossprod(col)))
    # Calculate of the reflection vector a-r
    a_r <- col; a_r[1] <- a_r[1] + sign(a_r[1]) * norm1
    # beta = 2 / ||a-r||^2  
    beta <- 2 / drop(crossprod(a_r))
    # update matrix Q (trailing matrix only) by Householder reflection
    Q[,k:n] <- Q[,k:n] - tcrossprod(Q[, k:n] %*% a_r, beta * a_r)
    # update matrix A (trailing matrix only) by Householder reflection
    A[k:n, k:p] <- A[k:n, k:p] - tcrossprod(beta * a_r, crossprod(A[k:n,k:p], a_r))
    }

  if (complete) {
     A[lower.tri(A)] <- 0
     return(list(Q = Q, R = A))
     }
  else {
    R <- A[1:p, ]; R[lower.tri(R)] <- 0
    return(list(Q = Q[,1:p], R = R))
    }
  }

Теперь давайте проведем тест:

X <- structure(c(0.8147, 0.9058, 0.127, 0.9134, 0.6324, 0.0975, 0.2785, 
0.5469, 0.9575, 0.9649, 0.1576, 0.9706, 0.9572, 0.4854, 0.8003
), .Dim = c(5L, 3L))

#       [,1]   [,2]   [,3]
#[1,] 0.8147 0.0975 0.1576
#[2,] 0.9058 0.2785 0.9706
#[3,] 0.1270 0.5469 0.9572
#[4,] 0.9134 0.9575 0.4854
#[5,] 0.6324 0.9649 0.8003

Сначала для версии с тонким QR:

## thin QR factorization
myqr(X)

#$Q
#            [,1]       [,2]        [,3]
#[1,] -0.49266686 -0.4806678  0.17795345
#[2,] -0.54775702 -0.3583492 -0.57774357
#[3,] -0.07679967  0.4754320 -0.63432053
#[4,] -0.55235290  0.3390549  0.48084552
#[5,] -0.38242607  0.5473120  0.03114461
#
#$R
#          [,1]       [,2]       [,3]
#[1,] -1.653653 -1.1404679 -1.2569776
#[2,]  0.000000  0.9660949  0.6341076
#[3,]  0.000000  0.0000000 -0.8815566

Теперь полная версия QR:

## full QR factorization
myqr(X, complete = TRUE)

#$Q
#            [,1]       [,2]        [,3]       [,4]       [,5]
#[1,] -0.49266686 -0.4806678  0.17795345 -0.6014653 -0.3644308
#[2,] -0.54775702 -0.3583492 -0.57774357  0.3760348  0.3104164
#[3,] -0.07679967  0.4754320 -0.63432053 -0.1497075 -0.5859107
#[4,] -0.55235290  0.3390549  0.48084552  0.5071050 -0.3026221
#[5,] -0.38242607  0.5473120  0.03114461 -0.4661217  0.5796209
#
#$R
#          [,1]       [,2]       [,3]
#[1,] -1.653653 -1.1404679 -1.2569776
#[2,]  0.000000  0.9660949  0.6341076
#[3,]  0.000000  0.0000000 -0.8815566
#[4,]  0.000000  0.0000000  0.0000000
#[5,]  0.000000  0.0000000  0.0000000

Теперь давайте проверим стандартный результат, возвращаемый qr.default:

QR <- qr.default(X)

## thin R factor
qr.R(QR)
#          [,1]       [,2]       [,3]
#[1,] -1.653653 -1.1404679 -1.2569776
#[2,]  0.000000  0.9660949  0.6341076
#[3,]  0.000000  0.0000000 -0.8815566

## thin Q factor
qr.Q(QR)
#            [,1]       [,2]        [,3]
#[1,] -0.49266686 -0.4806678  0.17795345
#[2,] -0.54775702 -0.3583492 -0.57774357
#[3,] -0.07679967  0.4754320 -0.63432053
#[4,] -0.55235290  0.3390549  0.48084552
#[5,] -0.38242607  0.5473120  0.03114461

## full Q factor
qr.Q(QR, complete = TRUE)
#            [,1]       [,2]        [,3]       [,4]       [,5]
#[1,] -0.49266686 -0.4806678  0.17795345 -0.6014653 -0.3644308
#[2,] -0.54775702 -0.3583492 -0.57774357  0.3760348  0.3104164
#[3,] -0.07679967  0.4754320 -0.63432053 -0.1497075 -0.5859107
#[4,] -0.55235290  0.3390549  0.48084552  0.5071050 -0.3026221
#[5,] -0.38242607  0.5473120  0.03114461 -0.4661217  0.5796209

Так что наши результаты верны!