Qr разложение матрицы: QR-разложение матрицы. Вычисление решения задачи наименьших… | by Iuliia Averianova | NOP::Nuances of Programming

QR. Расчет и использование QR-разложения | by Risto Hinno

Вычисление и использование QR-разложения

Введение

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

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

Разложение QR имеет следующую формулу:

A = QR, где:

• A исходная матрица, которую мы хотим разложить
• Q ортогональная матрица
• R основная треугольная матрица 90 цель довольно проста, разложить матрицу на матрицы Q и R. Чтобы найти ортогональную матрицу Q, мы могли бы использовать процесс Грама-Шмидта. Этот процесс берет входную матрицу и делает столбцы ортогональными (перпендикулярными друг другу). На следующем графике мы превращаем векторы a и b в ортогональные векторы v1 и v2.
  Ортогонализация векторов

  Для создания ортогональных векторов мы могли бы использовать следующую информацию. Вектор a является первым вектором, и мы можем сопоставить его с v1. v2 должен быть ортогонален вектору v1, поэтому мы должны спроецировать вектор b на v2, чтобы сделать b ортогоналом a (v1). Проекция вектора b на a (v1) определяется как:

  Информация, которую мы имеем, находится на следующем графике:

  Проецирование векторов

  Теперь нам нужно спроецировать b на v2, чтобы получить ортогональность между a и b. Для этого мы могли бы использовать e, которая является нашей ошибкой проекции и может быть рассчитана по следующей формуле:

  Следует отметить, что e имеет ту же длину и то же направление, что и v2. Мы можем заключить:

  Это в основном то, как работает процесс Грама-Шмидта. Для более чем двух измерений (предположим, что у нас есть векторы a, b и c) мы могли бы рассчитать процесс следующим образом (источник): v2 и т. д., если у нас больше измерений). Теперь мы знаем, как получить матрицу Q.

  Для вычисления матрицы R мы можем использовать следующие изменения в формуле QR:

  A = QR,

  QᵗA = QᵗQR (поскольку Q ортогонален Qᵗ=inverse(Q))

  QᵗA = R, таким образом, R = QᵗA

  In QR-разложение Python с использованием процесса Грама-Шмидта может быть реализовано следующим образом:

  Пример использования и вывода:

   A = np.array([[60, 91, 26], [60, 3, 75], [45, 90 , 31]], dtype='float') 
   Q, R = qr_factorization(A) 
   Q 
   #[[ 0,62469505, 0,71080736, -0,63928383], 
   # [0,62469505, 0,02343321, 0,75368929], 
   # [0,46852129, 0,70299629, -0,15254061]] 
   R 
   # [96,04686356, 77,6180, 
  , 
  , 
  , 
  , 0,6180.  0. , 0. , 35.17655816]]np.allclose(A, Q@R) 
   #True 

  Как обычно, это не единственный способ вычисления QR-разложения. Далее давайте посмотрим, как использовать отражения Хаусхолдера для QR-разложения.

  Размышления домохозяев

  Расчет Грама-Шмидта может быть численно неустойчивым, а это означает, что небольшие изменения на входе могут привести к относительно большим изменениям на выходе (источник). Более устойчивым способом является использование отражений Хаусхолдера. Проекты домохозяев векторируются через «зеркало». У нас есть вектор x, который мы хотим отразить в векторе Qx. Для отражения будем использовать ортогональную матрицу Q.

  Отражение Хаусхолдера

  Из предыдущего мы знаем, что проекция x на u равна:

  Из графика отражения Хаусхолдера мы можем видеть, что если мы дважды вычтем из x компонент, параллельный u, мы получим Qx (источник):

  Параллельный компонент оказывается нашей проекцией x на u, поэтому мы можем написать:

  Из-за ассоциативности матричного умножения:

  Мы можем написать:

  Что станет:

  И, наконец, мы имеем:

  Это способ получения матрицы Q. Полный процесс домохозяина состоит из следующих шагов:

  • выберите первый столбец (a) исходной матрицы A, которую мы хотим разложить:

   [[  12  , -51, 4], 
   [  6  , 167, -68], 
   [  - 4  , 24, -41]] 

  • вычислить вектор u, который перпендикулярен зеркалу, используя следующую формулу. Sign(a1) означает знак первого элемента в a и используется для предотвращения числовых ошибок отмены. e1 обозначает стандартный базисный вектор ([1, 0, 0, …], что означает, что все, кроме первой записи, станут равными нулю):
  Источник
  • получить матрицу Хаусхолдера H для столбца, используя предыдущую формулу:
  • получить Q, вычислив скалярное произведение Q = QH
  • обновить матрицу A с помощью матрицы домохозяев: A = HA i-й ряд и выполните предыдущие шаги. Пример выбора второго столбца обновленной матрицы A:

   [[-14., -21. ,  14. ] 
   [-0. ,  173,923,  -65,692] 
   [0.,  19.385,  -42.538]] 

  • Наконец, если процесс завершится, мы вернем Q и обновленную версию A в виде матриц Q и R.

  В python этот процесс реализован следующим образом:

  Если мы используем это с функцией, мы получаем:

   A = np.array([[60, 91, 26], [60, 3, 75], [45, 90, 31]], dtype='float') 
   Q, R = qr(a) 
   Q.round(3) 
   #[[-0,625, 0,355, -0,696], 
   # [-0,625, -0,762, 0,172], 
   # [-0,469, 0,542, 0,698]] 
   Р.раунд(3) 
  #[[-96.047,-100.888,-77.618],
  #[-0. , 78.813, -31.083], 
   # [ 0. , 0. , 16.469]]Q_np, R_np=np.linalg.qr(A) 
   np.allclose(Q, Q_np) 
   #True 
   np.allclose(R, R_np) 
   #True 

  Наша реализация дает такие же результаты, как и numpy.

  Теперь мы можем посмотреть, как QR-разложение можно использовать на практике.

  Решение линейных уравнений

  У нас есть матрицы A, x и B. Мы знаем A и B, но хотели бы найти x. Мы можем использовать некоторую алгебру (источник):

  • Ax = B
  • QRx = B
  • Rx = QᵗB (поскольку Qᵗ = inverse(Q), так как Q ортогонален)

  Решение последнего уравнения становится тривиальным с использованием обратной подстановки (см. дополнительную информацию здесь), потому что R – верхняя треугольная матрица. Простой пример:

   A = np.array([ 
   [2., 1., 1.], 
   [1., 3., 2.], 
   [1., 0., 0]]) 
   B = np.array([4., 5., 6.]) 
   #decompose 
   Q, R = np.linalg.qr(A) 
   #Q.TB 
   y = Q.T @ B 
   #решить x, code для back_substitution отсюда https://gist.github.com/RRisto/2cddc23d829877248915d66ecbe09a39 
   back_substitution(R, y) 
   #[ 6., 15., -23.]#using np Solr 
   np.linalg.solve(R, y) 
   #[ 6., 15., -23.] 

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

  Нахождение собственных значений и собственных векторов

  Нахождение собственных значений и векторов не очень сложно. Мы можем использовать следующую процедуру (источник):

  • создать матрицу X из матрицы A и разложить с помощью QR-разложения:

   Q, R = np.linalg.qr(X) 

  • вычисление произведения Q и предыдущих произведений (pQ инициализируется как диагональная матрица):

   pQ = pQ @ Q 

  • вычисление нового X с использованием R и Q:

   X = R @ Q 

  • повторить процесс n раз

  Собственные значения будут диагональными элементами X, а собственные векторы находятся в матрице pQ.

  Простая реализация Python для предыдущего алгоритма:

  Использование этой функции дает нам:

   A = np.array([[2., 1., 1.],[1., 3., 2.],[1., 0., 0]]) 
   find_eig_qr(A) 
   #eigenvalues 
   #[ 3.91222918,  1.28646207, -0.19869124] 
   #eigenvectors 
   #[[ 0.51233387, -0.72297297, -0.46349121], 
   # [ 0.84874276,  0.50856627,  0.144], 
   # [ 0.13095702, -0.46762212,  0.87417379]]#compare with реализация numpy 
   np. linalg.eig(A) 
   #eigenvalues ​​
   #[ 3.91222918, 1.28646207, -0.19869124] 
   #eigenvectors 
   #[[-0.51233387, -0.51118241, -0.118045], 
   # [-0.84874276, 0.76210326, -0.48349198], 
   # [-0.13095702, -0.39735522, 0.85856551]] 

  Мы видим, что собственные векторы отличаются от реализации, предоставленной numpy. Это очень простая реализация, но дает довольно хорошие результаты.

  Заключение

  QR-разложение может быть рассчитано с использованием процесса Грама-Шмидта и проекций Хаусхолдера, последние более стабильны в численном отношении. Вариантов этих алгоритмов больше, здесь я представил самые простые способы их использования. Разложение QR можно использовать для решения линейных уравнений и нахождения собственных значений и собственных векторов.

  Ссылки

  Поиск собственных векторов матриц с использованием QR-разложения, stats.stackexchange.com

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

  Линейные уравнения с Python: QR-разложение, b2bvoice. com

  0909 0

  QR-разложение, rosetta.com Преобразование Хаусхолдера в численной линейной алгебре, Джон Керл

  Почему Хаусхолдер в вычислительном отношении более устойчив, чем модифицированный Грам-Шмидт? math.stackexchange.com

  4. Разложение QR — количественная экономика с Python

  Количественная экономика с Python

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

  Томас Дж. Сарджент и Джон Стахурски

  4.1. Обзор

  В этой лекции описывается QR-разложение и его связь с

  • Ортогональная проекция и метод наименьших квадратов

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

  • Собственные значения и собственные векторы

  Мы напишем код на Python, чтобы закрепить наше понимание.

  4.2. Матричная факторизация

  QR-разложение (также называемое QR-факторизацией) матрицы — это разложение матрицы на произведение ортогональной матрицы и треугольной матрицы.

  QR-разложение вещественной матрицы \(A\) принимает вид

  \[ А=QR \]

  где

  Мы будем использовать процесс Грама-Шмидта для вычисления QR-разложения

  Поскольку это очень поучительно, мы напишем свой собственный код Python для выполнения этой работы

  4.

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

  Начнем с матрицы квадрата \(A\).

  Если квадратная матрица \(A\) невырожденна, то факторизация \(QR\) единственна.

  С прямоугольной матрицей \(A\) мы познакомимся позже.

  На самом деле наш алгоритм будет работать с прямоугольником \(A\), который не является квадратом.

  4.3.1. Процесс Грама-Шмидта для квадрата \(A\)

  Здесь мы применяем процесс Грама-Шмидта к столбцам матрицы \(A\).

  В частности, пусть

  \[ A = \left[ \begin{array}{c|c|c|c} a_1 & a_2 & \cdots & a_n \end{array} \right] \]

  Пусть \(|| · ||\) обозначает L2-норму.

  Алгоритм Грама-Шмидта многократно комбинирует следующие два шага в определенном порядке u_1=a_1, \ \ \ e_1=\frac{u_1}{||u_1||} \]

  Мы организуем сначала для вычисления \(u_2\), а затем нормализовать для создания \(e_2\):

  \[ u_2=a_2-(a_2· e_1)e_1, \ \ \ e_2=\frac{u_2}{||u_2||} \]

  Мы предлагаем читателю проверить, что \(e_1\) ортогонален \(e_2\), проверив, что \(e_1 \cdot e_2 = 0\).

  Процедура Грама-Шмидта продолжает повторяться.

  Таким образом, для \(k = 2, \ldots, n-1\) мы строим

  \[ u_{k+1}=a_{k+1}-(a_{k+1}· e_1)e_1-\cdots-(a_{k+1}· e_k)e_k, \ \ \ e_{k+1} =\frac{u_{k+1}}{||u_{k+1}||} \]

  Здесь \((a_j \cdot e_i)\) можно интерпретировать как линейный метод наименьших квадратов коэффициент регрессии \(a_j\) на \(e_i\)

  • это внутренний продукт \(a_j\) и \(e_i\), деленный на внутренний продукт \(e_i\), где \(e_i \cdot e_i = 1\), как заверила нас нормализация .

  • этот коэффициент регрессии интерпретируется как ковариация деленная на дисперсию

  Можно проверить, что

  \[\begin{split} A= \left[ \begin{array}{c|c|c|c} a_1 & a_2 & \cdots & a_n \end{array} \right]= \left[ \begin{array}{c|c|c|c} e_1 & e_2 & \cdots & e_n \end{массив} \right] \left[ \begin{matrix} a_1·e_1 & a_2·e_1 & \cdots & a_n·e_1\\ 0 & a_2·e_2 & \cdots & a_n·e_2 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_n·e_n \end{matrix} \right] \конец{разделить}\]

  Таким образом, мы построили разложение

  \[ А = Q р \]

  где

  \[ Q = \left[ \begin{array}{c|c|c|c} a_1 & a_2 & \cdots & a_n \end{array} \right]= \left[ \begin{array}{c|c|c|c} e_1 & e_2 & \cdots & e_n \end{массив} \right] \]

  и

  \[\begin{split} R = \left[ \begin{matrix} a_1·e_1 & a_2·e_1 & \cdots & a_n·e_1\\ 0 & a_2·e_2 & \cdots & a_n·e_2 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_n·e_n \end{matrix} \right] \конец{разделить}\]

  4.

  3.2. \(A\) не квадратная

  Теперь предположим, что \(A\) является матрицей \(n \times m\), где \(m > n\).

  Тогда \(QR\) разложение равно

  \[\begin{split} A = \left[ \begin{array}{c|c|c|c} a_1 & a_2 & \cdots & a_m \end{array} \right]=\left[ \begin{array}{c|c|c |c} e_1 & e_2 & \cdots & e_n \end{массив} \right] \left[ \begin{matrix} a_1·e_1 & a_2·e_1 & \cdots & a_n·e_1 & a_{n+1}\cdot e_1 & \cdots & a_{m}\cdot e_1 \\ 0 & a_2·e_2 & \cdots & a_n·e_2 & a_{n+1}\cdot e_2 & \cdots & a_{m}\cdot e_2 \\ \vdots & \vdots & \ddots & \quad \vdots & \ vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_n·e_n & a_{n+1}\cdot e_n & \cdots & a_{m}\cdot e_n \end{matrix} \right] \конец{разделить}\]

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

  \[\begin{align*} a_1 & = (a_1\cdot e_1) e_1 \cr a_2 & = (a_2\cdot e_1) e_1 + (a_2\cdot e_2) e_2 \cr \vdots & \quad \vdots \cr a_n & = (a_n\cdot e_1) e_1 + (a_n\cdot e_2) e_2 + \cdots + (a_n \cdot e_n) e_n \cr a_{n+1} & = (a_{n+1}\cdot e_1) e_1 + (a_{n+1}\cdot e_2) e_2 + \cdots + (a_{n+1}\cdot e_n) e_n \ кр \vdots & \quad \vdots \cr a_m & = (a_m\cdot e_1) e_1 + (a_m\cdot e_2) e_2 + \cdots + (a_m \cdot e_n) e_n \cr \end{выравнивание*}\]

  4.

  4. Some Code

  Теперь давайте напишем самодельный код Python для реализации декомпозиции QR путем развертывания описанного выше процесса Грама-Шмидта.

   импортировать numpy как np
  из scipy.linalg импортировать qr
   

   по определению QR_Decomposition(A):
      n, m = A.shape # получить форму A
      Q = np.empty((n, n)) # инициализируем матрицу Q
      u = np.empty((n, n)) # инициализируем матрицу u
      и[:, 0] = А[:, 0]
      Q[:, 0] = u[:, 0] / np.linalg.norm(u[:, 0])
      для i в диапазоне (1, n):
          и[:, я] = А[:, я]
          для j в диапазоне (i):
              u[:, i] -= (A[:, i] @ Q[:, j]) * Q[:, j] # получить каждый вектор u
          Q[:, i] = u[:, i] / np.linalg.norm(u[:, i]) # вычисляем каждый вектор
      R = np.zeros ((n, m))
      для я в диапазоне (n):
          для j в диапазоне (i, m):
              R[i, j] = A[:, j] @ Q[:, i]
      возврат Q, R
   

  Предыдущий код хорош, но может быть полезен при некоторых дополнительных действиях.

  Мы хотим сделать это, потому что позже в этой записной книжке мы хотим сравнить результаты использования нашего самодельного кода выше с кодом для QR, который предоставляет пакет Python scipy .

  Могут быть различия в знаках между матрицами \(Q\) и \(R\), созданными различными алгоритмами вычисления.

  Все они являются допустимыми QR-разложениями из-за того, что различия в знаках сокращаются при вычислении \(QR\).

  Однако, чтобы сделать результаты нашей самодельной функции и модуля QR в scipy сопоставимыми, давайте потребуем, чтобы \(Q\) имел положительные диагональные записи.

  Мы делаем это, корректируя знаки столбцов в \(Q\) и строк в \(R\) соответствующим образом.

  Для этого мы определим пару функций.

   определение diag_sign(A):
      «Вычислить знаки диагонали матрицы А»
      D = np.diag (np.sign (np.diag (A)))
      вернуть D
  определение настройки_знака (Q, R):
      """
      Отрегулируйте знаки столбцов в Q и строк в R, чтобы
      наложить положительную диагональ Q
      """
      D = диагностический_знак (Q)
      Q[:, :] = Q@D
      R[:, :] = D@R
      возврат Q, R
   

  4.5. Пример

  Теперь давайте сделаем пример.

   A = np.array([[1.0, 1.0, 0.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0]])
  # A = np.array([[1.0, 0.5, 0.2], [0.5, 0.5, 1.0], [0.0, 1.0, 1.0]])
  # A = np.array([[1.0, 0.5, 0.2], [0.5, 0.5, 1.0]])
  А
   

   массив([[1., 1., 0.],
         [1., 0., 1.],
         [0., 1., 1.]])
   

   Q, R = Adjust_sign(*QR_Decomposition(A))
   

  Массив

   ([[ 0,70710678, -0,40824829, -0,57735027],
         [0,70710678, 0,40824829, 0,57735027],
         [ 0. , -0,81649658, 0,57735027]])
   

  массив

  ([[ 1.41421356, 0.70710678, 0.70710678],
         [ 0. , -1,22474487, -0,40824829],
         [ 0. , 0. , 1.15470054]])
   

  Давайте сравним результаты с тем, что производит пакет scipy

   Q_scipy, R_scipy = Adjust_sign(*qr(A))
   

   print('Наш Q: \n', Q)
  печать('\п')
  print('Scipy Q:\n', Q_scipy)
   

   Наш вопрос:
   [[ 0,70710678 -0,40824829 -0,57735027]
   [0,70710678 0,40824829 0,57735027]
   [ 0. -0,81649658 0,57735027]]
  Сципи В:
   [[ 0,70710678 -0,40824829 -0,57735027]
   [0,70710678 0,40824829 0,57735027]
   [ 0.  -0,81649658 0,57735027]]
   

   print('Наш R: \n', R)
  печать('\п')
  print('Scipy R:\n', R_scipy)
   

   Наш Р:
   [[ 1,41421356 0,70710678 0,70710678]
   [ 0. -1,22474487 -0,40824829]
   [ 0. 0. 1.15470054]]
  Сципи Р:
   [[ 1,41421356 0,70710678 0,70710678]
   [ 0. -1,22474487 -0,40824829]
   [ 0. 0. 1.15470054]]
   

  Приведенные выше результаты дают нам хорошие новости: наша самодельная функция согласуется с тем, что scipy производит.

  Теперь давайте сделаем QR-разложение для прямоугольной матрицы \(A\), которая равна \(n \times m\) с \(т > п\).

   A = np.array([[1, 3, 4], [2, 0, 9]])
   

   Q, R = Adjust_sign(*QR_Decomposition(A))
  Q, Р
   

   (массив([[ 0,4472136 , -0,89442719],
          [0,89442719, 0,4472136]]),
   массив([[ 2.23606798, 1.34164079, 9.8386991],
          [ 0. , -2,68328157, 0,4472136 ]]))
   

   Q_scipy, R_scipy = Adjust_sign(*qr(A))
  Q_scipy, R_scipy
   

   (массив([[ 0,4472136 , -0,89442719],
          [0,89442719, 0,4472136]]),
   массив([[ 2. 23606798, 1.34164079, 9.8386991],
          [ 0. , -2,68328157, 0,4472136 ]]))
   

  4.6. Использование разложения QR для вычисления собственных значений

  Теперь немного полезного об алгоритме QR.

  Следующие итерации QR-разложения могут быть использованы для вычисления собственных значений из квадратных 9матрица 0004 \(A\).

  Вот алгоритм:

  1. Установить \(A_0 = A\) и сформировать \(A_0 = Q_0 R_0\)

  2. Форма \(A_1 = R_0 Q_0 \) . Обратите внимание, что \(A_1\) похож на \(A_0\) (легко проверить) и поэтому имеет те же собственные значения.

  3. Сформируйте \(A_1 = Q_1 R_1\) (т.е. сформируйте \(QR\) разложение \(A_1\)).

  4. Форма \(A_2 = R_1 Q_1 \) и затем \(A_2 = Q_2 R_2\) .

  5. Итерация до сходимости.

  6. Вычислите собственные значения \(A\) и сравните их с диагональными значениями предельного \(A_n\), полученного в результате этого процесса.

  Примечание: этот алгоритм близок к одному из самых эффективных способов вычисления собственных значений!

  Давайте напишем код на Python, чтобы опробовать алгоритм

   def QR_eigvals(A, tol=1e-12, maxiter=1000):
      «Найдите собственные значения A, используя QR-разложение».
      A_old = np.copy(A)
      A_new = np.copy(A)
      разница = np.inf
      я = 0
      в то время как (diff > tol) и (i < maxiter):
          A_старый[:, :] = A_новый
          Q, R = QR_Decomposition(A_old)
          A_new[:, :] = R @ Q
          diff = np.abs(A_new - A_old).max()
          я += 1
      eigvals = np.diag (A_new)
      вернуть эйвалы
   

  Теперь давайте попробуем код и сравним результаты с тем, что дает нам scipy.linalg.eigvals

  Вот

   # поэкспериментируйте с одной случайной матрицей A
  A = np.random.random ((3, 3))
   

   отсортировано (QR_eigvals (A))
   

   [-0,15277926583719334, 0,1820813029984996, 1,0321330823543302]
   

  Сравните с пакетом scipy .

   отсортировано (np.linalg.eigvals (A))
   

   [(0,014651018580652257-0,1293521629631194j),
   (0,014651018580652257+0,1293521629631194j),
   (1.032133082354331+0j)]
   

  4.7. \(QR\) и PCA

  Существуют интересные связи между \(QR\) разложением и анализом главных компонент (PCA).

  Вот некоторые.

  1. Пусть \(X'\) будет \(k \times n\) случайной матрицей, где \(j\)-й столбец является случайным числом из \({\mathcal N}(\mu, \Sigma)\), где \(\mu\) - это \(k \times 1\) вектор средних значений, а \(\Sigma\) - это \(k \times к\) ковариационная матрица. Мы хотим \(n > > k\) — это «эконометрический пример».

  2. Форма \(X' = Q R \), где \(Q \) равно \(k \times k\), а \(R\) равно \(k \times n\).

  3. Сформируем собственные значения \( R R'\), т. е. вычислим \(R R' = \tilde P \Lambda \tilde P' \).

  4. Сформируйте \(X' X = Q \tilde P \Lambda \tilde P' Q'\) и сравните его с собственным разложением \( X'X = P \hat \Lambda P'\).