Разное

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

8.3.4. QR-разложение MathCAD 12 руководство

  • Системы линейных уравнений
  • 8.1. Хорошо обусловленные системы с квадратной матрицей
    • 8.1.1. Вычислительный блок Given / Find
    • 8.1.2. Функция lsolve
  • 8.2. Произвольные системы линейных уравнений
    • 8.2.1. Переопределенные системы
    • 8.2.2. Недоопределенные системы
    • 8.2.3. Вырожденные и плохо обусловленные системы
  • 8.3. Матричные разложения
    • 8.3.1. СЛАУ с треугольной матрицей
    • 8.3.2. Разложение Холецкого
    • 8.3.3. LU-разложение
    • 8.3.4. QR-разложение
    • 8.3.5. SVD-(сингулярное) разложение
  • 8.4. Собственные векторы и собственные значения матриц

Среди матричных разложений особую роль играют ортогональные преобразования, обладающие свойством сохранения нормы вектора. Напомним из курса линейной алгебры, что матрица Q называется ортогональной, если Q

TQ=I, где I — единичная матрица. Свойство сохранения нормы вектора при ортогональных преобразованиях, т. е. |Qx|=|x|, дает рецепт поиска псевдорешения вырожденных СЛАУ, а именно, замену исходной задачи минимизации невязки с «плохой» матрицей |Аx-b| задачей |QT(Ax-b)|, в которой матрица уже будет «хорошей» благодаря специальному построению матрицы Q. Таким образом, ортогональные разложения используются при решении любых систем (в том числе с прямоугольной матрицей А, причем как переопределенных, так и недоопределенных).

Одним из важнейших вариантов ортогональных разложений некоторой матрицы А является QR-разложение вида A=QR, где Q — ортогональная матрица, а R — верхняя треугольная матрица:

  •  qr (A) — QR-разложение:
  •  А — вектор или матрица любого размера.

Результатом действия функции qr (А) является матрица, составленная из матриц Q и R соответственно (подобно рассмотренному в предыдущем разделе LU-разложению).

Выделить матрицы QR-разложения несложно при помощи встроенной функции выделения подматрицы submatrix (листинг 8.22).

Листинг 8.22. QR-разложение сингулярной матрицы

Если бы исходная СЛАУ Aх=b не была вырожденной, то можно было сразу записать: QTQ-Rx=QTb, откуда следует (благодаря ортогональности матрицы Q): Rx=QTb. Так как матрица к — треугольная, решение данной системы получается по формулам прямого хода. Если использовать нашу пользовательскую функцию решения треугольной системы (см. разд. 8.3. Т), то результат исходной СЛАУ запишется в виде одной строки кода: trg(R,QTb) (соответствующий пример решения невырожденной СЛАУ вы найдете на компакт-диске).

Обратимся теперь к проблеме решения СЛАУ с сингулярной квадратной NxN или прямоугольной NxM матрицей А. Если ее ранг равен k (он может быть меньше N и м, как в листинге 8.22, где N=M=S, a k=2), то получающаяся треугольная матрица R имеет следующую структуру:


где R1 — верхняя треугольная матрица, a R2 — просто некоторая матрица, а нули обозначают, в общем случае, нулевые матрицы соответствующих размеров.

Если система вырожденная, то она имеет бесконечное множество псевдорешений (векторов, минимизирующих норму невязки). При помощи QR-разложения можно сразу выписать одно из них (правда, не обладающее минимальной нормой). В нашем примере последняя строка матрицы R содержит одни нули, поэтому последняя компонента вектора псевдорешения х может быть абсолютно любой. Если положить х2=0, то остальные компоненты х определятся из треугольной СЛАУ R1x=QTb, как это проиллюстрировано листингом 8.23.

Листинг 8.23. Поиск одного из псевдорешений вырожденной СЛАУ посредством QR-разложения (продолжение листинга 8.22 )

Для того чтобы выбрать из всего множества псевдорешений (минимизирующих невязку исходной СЛАУ) нормальное псевдорешение (т. е. обладающее минимальной нормой), необходимо решить соответствующую задачу минимизации (см. разд. 8.2.3). Если построено QR-разложение, сделать это намного легче. Если произвольную компоненту решения обозначить переменной у: х

2=у, то она определится из соответствующей задачи оптимизации (листинг 8. 24, 1—3 строки), а остальные составляющие самого решения х — из треугольной СЛАУ R1x=QTb-R2y. В последней строке листинга выводится полученное нормальное псевдорешение, а также его норма и соответствующая норма невязки (которые полезно сравнить с результатом прошлого листинга). Вспомогательный рис. 8.13 помогает понять структуру минимизируемой функции из листинга 8.24.

Листинг 8.24. Нормальное псевдорешение вырожденной СЛАУ (продолжение листингов 8.22 и 8.23)

Рис. 8.13. Норма псевдорешения в зависимости от у (продолжение листинга 8.24)


Резюмируя практические аспекты применения QR-разложения, надо отметить, что алгоритмы решения СЛАУ на его основе практически одинаковы как для хорошо обусловленных, так и для сингулярных систем.

ПРИМЕЧАНИЕ

На прилагаемом к книге компакт-диске вы найдете дополнительные примеры построения QR-разложения как для квадратной, так и прямоугольной матрицы А.

Нравится

Твитнуть

28. Алгоритм qr-разложения.

Ортогональные матрицы и матрицы плоского вращения.

модулю 1.

Это следует из неравенства Коши-Шварца:

Сохранение углов между векторами следует из равенства:

QR-разложение может быть осуществлено методами вращения и отражения.

Рассмотрим вращение вектора на плоскости.

Матрица вращения задается в виде: , – угол вращения.

Свойство ортогональной матрицы – сохранение угла между векторами.

Видно, что матрица вращения – ортогональная матрица:

Если принять, что или , то .

Рассмотрим систему линейных алгебраических уравнений второго порядка:

Найдем матрицу Q такую, что

, где

Рассмотрим систему уравнений с матрицей .

Плоской матрицей вращения называется матрица, имеющая вид:

Можно подтвердить, что матица Q является ортогональной матрицей с определителем, равным 1.

Применение указанной матрицы к i-му столбцу матрицы A: , дает вектор , имеющий в j-ой позиции 0. (Верхний индекс обозначает номер вектора).

Применяя к исходной матрице указанные плоские матрицы вращения получим матрицу:

С помощью указанных матриц вращения все элементы матрицы R ниже главной диагонали становятся равными нулю.

Для исключения соответствующих элементов, коэффициенты c и s определяются выражениями:

Произведение ортогональных матриц является ортогональной матрицей.

Чтобы решить задачу, матрица A дополняется матрицей , матрица

A является произвольной.

Учитывая, что ортогональное преобразование вектора невязки:

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

Решение задачи наименьших квадратов при , сводится к задаче решения системы алгебраических уравнений с верхней треугольной матрицей:

Чтобы применить метод QR-разложения к решению задачи наименьших квадратов, нужно привести матрицу A к квадратной форме:

матрица B – произвольная.

Исходное уравнение:

Матрица является квадратной . К этой системе можно применить метод QR-разложения.

Применяя метод вращения, уравнение записывается: ,

размерность вектора , размерность вектора ,

– ортогональная матрица, – верхняя треугольная матрица.

Разобьем матрицу R на блоки:

Умножая матрицу R справа на можем записать:

невязка (ошибка)

От неизвестных параметров зависит только первое слагаемое нормы невязки.

Минимальное значение этого слагаемого, если матрица A максимальный размер, определяется из уравнения: .

Таким образом, задача наименьших квадратов решается в два этапа.

На первом этапе осуществляется QR-разложение расширенной матрицы и определяются ее подматрицы и .

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

Матрица Q является ортогональной матрицей, т. е. матрицей, транспонирование которой совпадает с обратной матрицей.

Матрица R – верхняя треугольная матрица, решение которой осуществляется методом обратной подстановки.

scipy.linalg.qr — Руководство по SciPy v1.10.0

scipy.linalg.qr(
a , overwrite_a=False , lwork=None , mode=’full’ , pivoting=False , check_finite=True 5 9000) [016]

Вычисление QR-разложения матрицы.

Вычислить разложение A = Q R где Q унитарно/ортогонально и R верхний треугольный.

Параметры:
а (M, N) array_like

Матрица для разложения

overwrite_a bool, необязательно

Переписываются ли данные в a (может улучшить производительность, если overwrite_a устанавливается в True путем повторного использования существующих входных данных структуру, а не создавать новую. )

lwork int, необязательный

Размер рабочего массива, lwork >= a.shape[1]. Если None или -1, оптимальный размер вычисляется.

режим {‘полный’, ‘r’, ‘экономичный’, ‘необработанный’}, необязательный

Определяет, какая информация должна быть возвращена: либо Q, либо R («полный», по умолчанию), только R («r») или оба Q и R, но вычисляются в эконом-размер («экономичный», см. Примечания). Окончательный вариант «сырой» (добавлено в SciPy 0.11) заставляет функцию возвращать две матрицы (Q, TAU) во внутреннем формате, используемом LAPACK.

поворот bool, необязательно

Должна ли факторизация включать поворот для раскрытия ранга qr разложение. При повороте вычислить разложение A P = Q R , как указано выше, но где P выбран так, что диагональ R не возрастает.

check_finite bool, необязательный

Проверять, содержит ли входная матрица только конечные числа. Отключение может дать прирост производительности, но может привести к проблемам (сбои, незавершение), если входные данные содержат бесконечности или NaN.

Возвращает:
Q число с плавающей запятой или комплексный ndarray

Форма (M, M) или (M, K) для mode='economic' . Не возвращено если режим='r' .

R с плавающей запятой или комплексный ndarray

Форма (M, N) или (K, N) для mode='economic' . К = мин(М, Н) .

P int ndarray

Формы (N,) для поворота = True . Не возвращается, если поворот=ложь .

Поднимает:
LinAlgError

Возникает при сбое декомпозиции

Примечания

Это интерфейс для подпрограмм LAPACK dgeqrf, zgeqrf, dorgqr, zungqr, dgeqp3 и zgeqp3.

Если режим = экономичный , формы Q и R вместо этого будут (M, K) и (K, N) из (M,M) и (M,N), где K=min(M,N) .

Примеры

 >>> импортировать numpy как np
>>> из scipy import linalg
>>> rng = np.random.default_rng()
>>> a = rng.standard_normal((9, 6))
 
 >>> q, r = linalg.qr(a)
>>> np.allclose(a, np.dot(q, r))
Истинный
>>> q.shape, r.shape
((9, 9), (9, 6))
 
 >>> r2 = linalg.qr(a, mode='r')
>>> np.allclose(r, r2)
Истинный
 
 >>> q3, r3 = linalg.qr(a, режим='экономический')
>>> q3.форма, r3.форма
((9, 6), (6, 6))
 
 >>> q4, r4, p4 = linalg.qr(a, pivoting=True)
>>> d = np.abs(np.diag(r4))
>>> np.all(d[1:] <= d[:-1])
Истинный
>>> np.allclose(a[:, p4], np.dot(q4, r4))
Истинный
>>> q4.форма, r4.форма, p4.форма
((9, 9), (9, 6), (6,))
 
 >>> q5, r5, p5 = linalg.qr(a, mode='economic', pivoting=True)
>>> q5.форма, r5.форма, p5.форма
((9, 6), (6, 6), (6,))

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

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

QR-факторизация (расширенный)

Мы видели, что для вычисления LSE необходимо инвертировать матрицу. В предыдущих разделах мы использовали функцию решить . Однако решение не является стабильным решением. При кодировании вычислений LSE мы используем разложение QR. 93) colnames(X) <- c("Перехват","x","x2","x3") бета <- матрица (с (1,1,1,1),4,1) set.seed(1) y <- X%*%beta+rnorm(n,sd=1)

Стандартная функция R для инверсии выдает ошибку:

solve(crossprod(X)) %*% crossprod(X,y)

Чтобы понять, почему это происходит, посмотрите

 варианты (цифры = 4)
log10 (кросспродукт (X))
 
 ## Перехват x x2 x3
## Пересечение 1,699 4,098 6,625 9,203
## x 4,098 6,625 9,203 11,810
## x2 6,625 9.203 11.810 14.434
## x3 9,203 11,810 14,434 17,070

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

Факторизация

QR-факторизация основана на математическом результате, который говорит нам, что мы можем разложить любую матрицу полного ранга следующим образом:

с:

  • матрица с
  • верхняя треугольная матрица.

Верхнетреугольные матрицы очень удобны для решения системы уравнений.

Пример верхней треугольной матрицы

В приведенном ниже примере матрица слева является верхнетреугольной: она имеет только 0 под диагональю. Это значительно облегчает решение системы уравнений:

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

Нахождение LSE с помощью QR

Если мы перепишем уравнения LSE, используя вместо у нас есть:

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

Теперь мы готовы найти LSE, используя разложение QR. Чтобы решить:

Мы используем обратное решение , которое использует верхнетреугольную природу .

 QR <- QR(X)
Q <- qr.Q( QR )
R <- qr.R( QR )
(бетахат <- обратное решение (R, crossprod (Q, y)) )
 
 ## [1]
## [1,] 0,9038
## [2,] 1,0066
## [3,] 1.0000
## [4,] 1.0000

На практике ничего из этого делать не нужно благодаря встроенной функции solv.qr :

 QR <- qr(X)
(бетахат <- решить.qr(QR, y))
 
 ## [1]
## Перехват 0,9038
## х 1,0066
## x2 1.0000
## x3 1.0000
Подобранные значения

Эта факторизация также упрощает вычисление подходящих значений:

В R мы просто делаем следующее: 92)/дф varbeta <- sigma2*chol2inv(qr.R(QR)) SE <- sqrt(diag(varbeta)) cbind (бетахат, SE) 

 ## ЮВ
## Перехват 0,9038 4.508e-01
## x 1,0066 7,858e-03
## x2 1.0000 3.662e-05
## x3 1.0000 4.802e-08

Это дает нам те же результаты, что и функция lm .

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *