Наука о данных и разложение матриц 

Специалистам по данным стоит хорошо знать несколько разложений матриц, потому что они помогают находить методы для актуальных вычислений и оценки результатов для моделей и алгоритмов, с которыми мы работаем. В некоторых случаях конкретной формой разложения является алгоритм (например, в методе главных компонент и сингулярном разложении). И во всех случаях разложение матриц помогает развивать интуицию и способность к анализу.

QR—разложение — одно из самых полезных. Ему находится множество важных применений в науке о данных, статистике и анализе данных. Одно из подобных применений — вычисление решения задачи наименьших квадратов. 

Содержание статьи

  • Повторение проблемы наименьших квадратов.
  • Описание с QR-разложения.
  • Решение задачи наименьших квадратов с помощью QR-разложения.
  • Реализация вычисления QR-разложения на R и C++ и сравнение реализаций.

Задача наименьших квадратов

QR-разложение позволяет вычислить решение задачи наименьших квадратов. Подчеркиваю, именно вычислить, поскольку обычный метод наименьших квадратов дает нам закрытое решение в форме нормальных уравнений. Это замечательно, но, если нам нужно найти актуальное числовое решение, этот метод не подойдет.

Вспомним задачу наименьших квадратов. Нужно решить уравнение ниже: 

Проблема состоит в том, что не существует решения для β, потому что обычно, если у нас больше наблюдений, чем переменных, X не имеет обратного значения, следовательно, вычисление ниже невозможно: 

Вместо этого попробуем найти некоторое β̂, решающее уравнение неидеально, но с минимально возможной ошибкой. Один из способов— минимизировать следующую целевую функцию, являющуюся функцией от β̂.

Минимизация этой суммы квадратов отклонений и дает имя задаче наименьших квадратов. Взятие производных по β̂ и приравнивание к нулю приведут к нормальным уравнениям и предоставят решение в замкнутой форме. 

Это один из способов. Но можно использовать линейную алгебру. И вот здесь QR-разложение подойдет как нельзя лучше. 

QR-разложение

Для начала давайте опишем это разложение. QR-разложение позволяет отобразить матрицу как произведение двух отдельных матриц Q и R.

Q  —  ортогональная матрица, а R  —  треугольная матрица.

Это означает, что:

Так как R квадратная, до тех пор, пока диагональные элементы не равны нулю, она также обратима. Если столбцы X линейно независимы, это всегда будет верным. Хотя, если в данных есть коллинеарность, все же будут возникать проблемы. Тем не менее — и в этом суть QR-разложения — прямоугольная и необратимая X может быть выражена как две обратимые матрицы! И вот это уже имеет смысл. 

Решение задачи наименьших квадратов с помощью QR-разложения.

Теперь, зная, что из себя представляет QR-разложение, решим задачу наименьших квадратов следующим образом: 

Следовательно:

Это означает, что все, что нужно сделать — это найти матрицу, обратную R, транспонировать Q и вычислить их произведение. Мы получим коэффициенты обычного метода наименьших квадратов. Нам даже не нужно вычислять ковариационную матрицу и обратную ей, что происходит в решении обычным методом наименьших квадратов. 

Реализация QR-разложения

Чтобы найти QR-коэффициенты матрицы, сначала найдем Q, используя процесс Грама-Шмидта, затем просто умножим исходную матрицу на транспонированную Q, чтобы найти R. Далее применим QR-разложение, используя функции, внедренные в R и C++. Позднее мы заглянем внутрь этих функций, чтобы рассмотреть весь процесс в деталях.

Вычисление коэффициентов

Я использую две функции: myQRR и myQRCpp, применяющие процесс Грама-Шмидта для QR-разложения. Одна функция написана на R, а другая на C++, она загружается в среду R через Rcpp. В конце статьи я сравню их производительность. 

library(Rcpp)
library(tidyverse)
library(microbenchmark)

sourceCpp("../source-code/myQRC.cpp")
source("../source-code/myQRR.R")

Начнем с маленького примера, в котором смоделируем y и X, а затем решим уравнение, используя QR-разложение. Также мы сможем провести двойную проверку QR-разложения — работает ли оно и возвращает ли смоделированный X. Вот смоделированные переменные ответа: 

y = rnorm(6)
y

## [1] 0.6914727 2.4810138 0.4049580 0.3117301 0.6084374 1.4778950

Вот данные, которые мы используем для определения коэффициентов наименьших квадратов. В нашем распоряжении 3 переменные: 

X = matrix(c(3, 2, 3, 2, -1, 4,
             5, 1, -5, 4, 2, 1,
             9, -3, 2 , -1, 8, 1), ncol = 3)
X

##      [,1] [,2] [,3]
## [1,]    3    5    9
## [2,]    2    1   -3
## [3,]    3   -5    2
## [4,]    2    4   -1
## [5,]   -1    2    8
## [6,]    4    1    1

Теперь, чтобы найти Q и R, используем myQRCpp.

  1. Видим, что R действительно верхнетреугольная матрица. 
Q = myQRCpp(X)$Q
R = t(Q) %*% X %>% round(14)

R

##          [,1]     [,2]      [,3]
## [1,] 6.557439 1.829983  3.202470
## [2,] 0.000000 8.285600  4.723802
## [3,] 0.000000 0.000000 11.288484

Q

##            [,1]        [,2]        [,3]
## [1,]  0.4574957  0.50241272  0.45724344
## [2,]  0.3049971  0.05332872 -0.37459932
## [3,]  0.4574957 -0.70450052  0.34218986
## [4,]  0.3049971  0.41540270 -0.34894183
## [5,] -0.1524986  0.27506395  0.63684585
## [6,]  0.6099943 -0.01403387 -0.07859294

2. Убеждаемся в ортогональности Q.

t(Q)%*%Q %>% round(14)

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

3. И в том, что QR действительно возвращает исходную матрицу X.

Q %*% R

##      [,1] [,2] [,3]
## [1,]    3    5    9
## [2,]    2    1   -3
## [3,]    3   -5    2
## [4,]    2    4   -1
## [5,]   -1    2    8
## [6,]    4    1    1

Теперь вычисляем актуальные коэффициенты.

beta_qr = solve(R) %*% t(Q) %*% y

beta_qr

##             [,1]
## [1,]  0.32297414
## [2,]  0.07255123
## [3,] -0.02764668

Чтобы проверить, верен ли ответ, сравним вычисленное значение β̂ с тем, что выдает функция lm.

coef(lm(y ~ -1 + ., data = data.frame(cbind(y,X))))

##          V2          V3          V4 
##  0.32297414  0.07255123 -0.02764668

Мы получили в точности то же решение, что и для рассчитанных коэффициентов. 

Реализация QR-разложения

Процесс Грама-Шмидта

Процесс Грама-Шмидта — это метод вычисления ортогональной матрицы Q, которая состоит из ортогональных или независимых единичных векторов и занимает такое же пространство, что и матрица X.

  • В качестве начального шага алгоритм включает в себя выбор вектора столбца X, скажем, x1 = u1.
  • Затем находим вектор, ортогональный u1, проецируя на него следующий столбец X, скажем, x2, и вычитая из него проекцию u2 = x2 − proj u1x2. Теперь у нас есть набор из двух ортогональных векторов. 
  • Следующий шаг — проделать то же самое, но вычитая сумму проекций на каждый вектор из набора ортогональных векторов uk.

Выразим это следующим образом:

Получив полный набор ортогональных векторов, просто делим каждый вектор на его нормаль и помещаем их в матрицу:

Зная Q, легко вычисляем R:

Реализация в R и C++

Конечно, в R существуют встроенные функции, осуществляющие QR-разложение. Так как алгоритм Грама-Шмидта является итеративным по своей природе, я решил реализовать его на C++ , являющемся хорошим инструментом для подобных задач, и сравнить его с подобной функцией R. Вот как выглядит моя версия R:

myQR = function(A){
  dimU = dim(A)
  U = matrix(nrow = dimU[1], ncol = dimU[2])
  U[,1] = A[,1]
  for(k in 2:dimU[2]){
    subt = 0
    j = 1
    while(j < k){
      subt = subt + proj(U[,j], A[,k])
      j = j + 1
    }
    U[,k] = A[,k] - subt
  }
  Q = apply(U, 2, function(x) x/sqrt(sum(x^2)))
  R = round(t(Q) %*% A, 10)
  return(list(Q = Q, R = R, U = U))
}

Очень точно. Внутри цикла for есть цикл while, и вызываемая проекционная функция также является функцией, написанной на R.

А вот как выглядит моя версия C++. Логика по сути та же, за исключением того, что есть другой цикл for, нормализующий ортогональные столбцы. 

// [[Rcpp::export]]
List myQRCpp(NumericMatrix A) {
  int a = A.rows();
  int b = A.cols();
  NumericMatrix U(a,b);
  NumericMatrix Q(a,b);
  NumericMatrix R(a,b);
  NumericMatrix::Column Ucol1 = U(_ , 0);
  NumericMatrix::Column Acol1 = A(_ , 0);
  
  Ucol1 = Acol1;
  
  for(int i = 1; i < b; i++){
    NumericMatrix::Column Ucol = U(_ , i);
    NumericMatrix::Column Acol = A(_ , i);
    NumericVector subt(a);
    int j = 0;
    while(j < i){
      NumericVector uj = U(_ , j);
      NumericVector ai = A(_ , i);
      subt = subt + projC(uj, ai);
      j++;
    }
    Ucol = Acol - subt;
  }
  for(int i = 0; i < b; i++){
    NumericMatrix::Column ui = U(_ , i);
    NumericMatrix::Column qi = Q(_ , i);
    
    double sum2_ui = 0;
    for(int j = 0; j < a; j++){
      sum2_ui = sum2_ui + ui[j]*ui[j];
    }
    qi = ui/sqrt(sum2_ui);
  }
  
  List L = List::create(Named("Q") = Q , _["U"] = U);
  return L;
}

Сравнение реализаций на R и C++ 

В дополнение к двум функциям выше у меня есть третья, идентичная R, только она называется projC, а не proj. Я назвал ее myQRC (projC написана на C++, а proj — на R). В противном случае у нас есть чистая C++ функция myQRCpp и чистая R функция myQR.

Чтобы сравнить, как быстро эти три функции выполняют QR-разложение, я поместил их в функцию QR_comp, которая вызывает функции и измеряет время выполнения каждой с одним и тем же аргументом матрицы. 

QR_comp = function(A){
  t0 = Sys.time()
  myQR(A)
  tQR = Sys.time() - t0
  t0 = Sys.time()
  myQRC(A)
  tQRC = Sys.time() - t0
  t0 = Sys.time()
  myQRCpp(A)
  tQRCpp = Sys.time() - t0
  
  return(data.frame(tQR = as.numeric(tQR), 
                    tQRC = as.numeric(tQRC),
                    tQRCpp = as.numeric(tQRCpp)))
}

Их производительность можно сравнить по сетке из случайных матриц n на m. Эти матрицы генерируются при вызове функции QR_comp.

grid = expand.grid(n = seq(10, 3010, 500), 
                   m = seq(50, 600, 50))

tvec = map2(grid$n, 
            grid$m, 
            ~QR_comp(matrix(runif(.x*.y), ncol = .y)))

Наконец, визуализируем.

plotly::ggplotly(
bind_rows(tvec) %>%
  gather("func","time") %>%
  mutate(n = rep(grid$n, 3),
         m = rep(grid$m, 3)) %>%
  ggplot(aes(m, n, fill = time)) + 
  geom_tile() + 
  facet_grid(.~func) +
  scale_fill_gradientn(colours = rainbow(9)) +
  theme(panel.background = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.y = element_text(angle = 35, size = 5),
        axis.text.x = element_text(angle = 30, size = 5)), width = 550, heigh = 400)

Очевидно, что чем больше задействован C++ , тем быстрее вычисление QR-разложения. Функция C++ решает его менее, чем за минуту, для матриц размером до 250 столбцов на 3000 строк или 600 столбцов на 500. R-функция в 2–3 раза медленнее.

Заключение

QR — это просто разложение матрицы, а метод наименьших квадратов просто одно из применений QR. Надеюсь, обсуждение выше демонстрирует, насколько важна и полезна линейная алгебра для науки о данных.

Читайте также:

Перевод статьи Ben Denis Shaffer: QR Matrix Factorization

Предыдущая статьяGo: как циклы преобразуются в ассемблерную программу?
Следующая статья64 совета на основе 50 лет опыта в разработке ПО