Численные методы линейной алгебры
Численные методы линейной алгебры (прикладная линейная алгебра) — это методы приближенного решения задач из области вычислительной математики и линейной алгебры. Целью дисциплины является разработка и анализ алгоритмов для численного решения матричных задач. Наиболее важными задачами являются решение систем линейных алгебраических уравнений и вычисление собственных значений.
Численная линейная алгебра исследует то, как можно использовать матричные операции для создания компьютерных алгоритмов, которые эффективно дают приблизительные ответы на задачи непрерывной математики. Численная линейная алгебра является видом линейной алгебры и представляет собой область численного анализа. Компьютеры используют арифметику с плавающей запятой и не могут точно представлять иррациональные данные, поэтому, когда компьютерный алгоритм применяется к матрице данных, он иногда может увеличить разность[англ.] между числом, хранящимся в компьютере, и истинным числом, приближением которого оно является. Численная линейная алгебра использует свойства векторов и матриц для разработки эффективных алгоритмов, которые минимизируют ошибку, вносимую компьютером.
Численная линейная алгебра направлена на решение задач непрерывной математики с использованием компьютеров конечной точности, поэтому ее приложения в естественных и общественных науках столь же обширны, как и приложения непрерывной математики. Она часто является фундаментальной частью инженерных и вычислительных проблем, таких как обработка изображений и сигналов, электросвязь, финансовые вычисления[англ.], материаловедение, структурная биология, интеллектуальный анализ данных, биоинформатика и гидродинамика. Особенно широко используются матричные методы в методах конечных разностей, методах конечных элементов и моделировании дифференциальных уравнений. Отмечая широкое применение численной линейной алгебры, Ллойд Н. Трефетен[англ.] и Дэвид Бау, III утверждают, что она «столь же фундаментальна для математических наук, как исчисление и дифференциальные уравнения»[1]:x, хотя это сравнительно небольшая область[2]. Поскольку многие свойства матриц и векторов также применимы к функциям и операторам, числовую линейную алгебру также можно рассматривать как тип функционального анализа, в котором особое внимание уделяется практическим алгоритмам[1]:ix.
Основные задачи в численной линейной алгебре включают получение матричных разложений, таких как сингулярное разложение, QR-разложение, LU-разложение или собственное разложение, которые затем можно использовать для решения общих задач линейной алгебры, таких как решение систем линейных уравнений, определение местоположения собственных значений или оптимизация методом наименьших квадратов. Основная задача численной линейной алгебры — разработка алгоритмов, которые не вносят ошибок при применении к реальным данным на компьютере с конечной точностью, часто достигается с помощью итерационных методов[англ.].
История
Численная линейная алгебра была разработана Джоном фон Нейманом, Аланом Тьюрингом, Джеймсом Х. Уилкинсоном, Алстоном Скоттом Хаусхолдером[англ.], Джорджем Форсайтом[англ.]* и Хайнцем Рутисхаузером[англ.], для применения самых ранних компьютеров к задачам непрерывной математики, таким как баллистические задачи и задачи на решение системы дифференциальных уравнений в частных производных[2]. Первая серьезная попытка свести к минимуму вычислительную ошибку при применении алгоритмов к реальным данным была сделана Джоном фон Нейманом и Германом Голдстайном в 1947 году[3]. Эта область расширялась по мере того, как технологии всё чаще позволяли исследователям решать сложные задачи на чрезвычайно больших высокоточных матрицах, а некоторые численные алгоритмы приобретали известность по мере того, как такие технологии, как параллельные вычисления, позволили применить практический подход к научным проблемам[2].
Разложения матриц
Блочная матрица
Для многих задач прикладной линейной алгебры полезно рассматривать матрицу как объединённый набор векторов-столбцов. Например, при решении линейной системы [math]\displaystyle{ x = A^{-1}b }[/math], вместо нахождения [math]\displaystyle{ x }[/math] как умножения [math]\displaystyle{ A^{-1} }[/math] и [math]\displaystyle{ b }[/math], лучше представить [math]\displaystyle{ x }[/math] как вектор коэффициентов в линейном разложении [math]\displaystyle{ b }[/math] в базисе, образованном столбцами [math]\displaystyle{ A }[/math][1]:8. Представление о матрицах как об объединённом наборе столбцов практично для целей матричных алгоритмов, в связи с тем, что матричные алгоритмы часто содержат два вложенных цикла: один по столбцам матрицы [math]\displaystyle{ A }[/math], а другой по строкам матрицы [math]\displaystyle{ A }[/math]. Например, для матриц [math]\displaystyle{ A^{m \times n} }[/math] и векторов [math]\displaystyle{ x^{n \times 1} }[/math] и [math]\displaystyle{ y^{m \times 1} }[/math], можно использовать разделение столбцов для вычисления [math]\displaystyle{ Ax+y }[/math] как
for j = 1:n
for i = 1:m
y(i) = A(i,j)x(j) + y(i)
end
end
Сингулярное разложение
Сингулярное разложение матрицы [math]\displaystyle{ A^{m \times n} = U \Sigma V^\ast }[/math], где [math]\displaystyle{ U }[/math] и [math]\displaystyle{ V }[/math] унитарные матрицы, а [math]\displaystyle{ \Sigma }[/math] диагональная. Элементы диагональной матрицы [math]\displaystyle{ \Sigma }[/math] называются сингулярными числами матрицы [math]\displaystyle{ A }[/math]. Так как сингулярные числа являются квадратными корнями из собственных значений матрицы [math]\displaystyle{ AA^\ast }[/math], существует тесная связь между разложением по сингулярным значениям и разложениями по собственным значениям. Это означает, что большинство методов вычисления разложения по сингулярным числам аналогичны методам собственных значений[1]:36; возможно, самый распространенный метод включает преобразование Хаусхолдера[1]:253.
QR-разложение
QR-разложение матрицы [math]\displaystyle{ A^{m \times n} }[/math] представляется произведением матрицы [math]\displaystyle{ Q^{m \times m} }[/math] на [math]\displaystyle{ R^{m \times n} }[/math] таким, что [math]\displaystyle{ A = QR }[/math], где [math]\displaystyle{ Q }[/math] — ортогональная матрица и [math]\displaystyle{ R }[/math] — верхняя треугольная матрица[1]:50,[4]:223. Два основных алгоритма для вычисления QR-разложения: процесс Грама-Шмидта и преобразование Хаусхолдера. QR-разложение часто используется для задач, связанных с методом наименьших квадратов, и задач на собственные значения (с помощью итеративного QR-алгоритма).
LU-разложение
LU-разложение матрицы [math]\displaystyle{ A }[/math] состоит из нижней треугольной матрицы [math]\displaystyle{ L }[/math] и верхней треугольной матрицы [math]\displaystyle{ U }[/math] такими, что [math]\displaystyle{ A = LU }[/math]. Матрица [math]\displaystyle{ U }[/math] находится с помощью метода Гаусса, которая включает в себя умножения [math]\displaystyle{ A }[/math] слева на набор матриц [math]\displaystyle{ M_1,\ldots,M_{n-1} }[/math] для формирования [math]\displaystyle{ M_{n-1} \cdots M_1 A = U }[/math], откуда [math]\displaystyle{ L = M_1^{-1} \cdots M_{n-1}^{-1} }[/math][1]:147[4]:96.
Спектральное разложение
Спектральное разложение матрицы [math]\displaystyle{ A^{m \times m} = X \Lambda X^{-1} }[/math], где столбцы [math]\displaystyle{ X }[/math] — собственные векторы матрицы [math]\displaystyle{ A }[/math], и [math]\displaystyle{ \Lambda }[/math] представляет собой диагональную матрицу, диагональные элементы которой — соответствующие собственные значения матрицы [math]\displaystyle{ A }[/math][1]:33. Не существует прямого метода нахождения спектрального разложения произвольной матрицы: поскольку невозможно написать программу, которая находит точные корни произвольного многочлена за конечное время, любой алгоритм нахождения собственных значений обязательно должен быть итерационным[1]:192.
Алгоритмы
Метод Гаусса
С точки зрения числовой линейной алгебры, метод Гаусса — это процедура разложения матрицы [math]\displaystyle{ A }[/math] в LU-разложение, которую метод Гаусса выполняет путем умножения [math]\displaystyle{ A }[/math] слева на последовательность матриц [math]\displaystyle{ L_{m-1} \cdots L_2 L_1 A = U }[/math] до тех пор, пока [math]\displaystyle{ U }[/math] не станет верхней треугольной, а [math]\displaystyle{ L }[/math] не станет нижней треугольной, где [math]\displaystyle{ L \equiv L_1^{-1}L_2^{-1} \cdots L_{m-1}^{-1} }[/math][1]:148. Наивные программы для метода Гаусса, как известно, очень вычислительно неустойчивы и дают огромные ошибки при применении к матрицам со многими значащими цифрами[2]. Самое простое решение — ввести опорный элемент[англ.], что даёт модифицированный стабильный алгоритм метода Гаусса[1]:151.
Решения линейных систем
Численная линейная алгебра обычно рассматривает матрицы как объединённый набор векторов-столбцов. Традиционный подход состоит в том, чтобы решить линейную систему [math]\displaystyle{ x = A^{-1}b }[/math], чтобы получить [math]\displaystyle{ x }[/math] как произведение [math]\displaystyle{ A^{-1} }[/math] и [math]\displaystyle{ b }[/math]. Вместо этого численная линейная алгебра интерпретирует [math]\displaystyle{ x }[/math] как вектор коэффициентов линейного разложения [math]\displaystyle{ b }[/math] в базисе, образованном столбцами [math]\displaystyle{ A }[/math][1]:8.
Для решения систем линейных уравнений могут применяться матричные уравнения. В зависимости от характеристик матрицы [math]\displaystyle{ A }[/math] и векторов [math]\displaystyle{ x }[/math] и [math]\displaystyle{ b }[/math], одни разложения могут быть проще, чем другие. можно использовать множество различных разложений, в зависимости от характеристик матрицы [math]\displaystyle{ A }[/math] и векторов [math]\displaystyle{ x }[/math] и [math]\displaystyle{ b }[/math], что может значительно упростить получение одного разложения по сравнению с другими. Если [math]\displaystyle{ A=QR }[/math] — QR-разложение матрицы [math]\displaystyle{ A }[/math], то [math]\displaystyle{ Rx = Q^\ast b }[/math][1]:54. Если [math]\displaystyle{ A = X\Lambda X^{-1} }[/math] — спектральное разложение матрицы [math]\displaystyle{ A }[/math], то мы стремимся найти [math]\displaystyle{ b }[/math] такой, что [math]\displaystyle{ b=Ax }[/math], с [math]\displaystyle{ b^\prime = X^{-1}b }[/math] и [math]\displaystyle{ x^\prime = X^{-1}x }[/math]. Откуда получаем [math]\displaystyle{ b^\prime = \Lambda x^\prime }[/math][1]:33. Наконец, если [math]\displaystyle{ A=LU }[/math] — LU-разложение матрицы [math]\displaystyle{ A }[/math], то [math]\displaystyle{ Ax=b }[/math] можно решить с помощью треугольных матриц [math]\displaystyle{ Ly=b }[/math] и [math]\displaystyle{ Ux = y }[/math][1]:147[4]:99.
Оптимизация методом наименьших квадратов
Подход, связанный с матричными разложениями, предлагает ряд способов решения линейной системы [math]\displaystyle{ r=b-Ax }[/math], где мы стремимся минимизировать [math]\displaystyle{ r }[/math], как в модели линейной регрессии. QR-алгоритм решает эту проблему, сначала определяя [math]\displaystyle{ y= Ax }[/math], а затем вычисляя тонкое QR-разложенние [math]\displaystyle{ A }[/math] и переходит к уравнению [math]\displaystyle{ \widehat{R}x = \widehat{Q}^\ast b }[/math]. Эта верхнетреугольная система может быть решена относительно [math]\displaystyle{ b }[/math]. Аналогичный алгоритм для решения МНК можно получить при помощи SVD-разложения. Вычисляя тонкое SVD-разложение [math]\displaystyle{ A = \widehat{U}\widehat{\Sigma}V^\ast }[/math], а затем вычисляя вектор [math]\displaystyle{ \widehat{U}^\ast b }[/math], мы сводим задачу наименьших квадратов к простой диагональной системе[1]:84. Тот факт, что решения методом наименьших квадратов могут быть получены с помощью QR-разложения и SVD-разложения, означает, что в дополнение к классическому методу нормальных уравнений[англ.] решения задач наименьших квадратов, эти задачи также могут быть решены методами, включающими алгоритм Грама-Шмидта и методы Хаусхолдера.
Обусловленность и вычислительная устойчивость
Рассмотрим функцию [math]\displaystyle{ f: X \to Y }[/math], где [math]\displaystyle{ X }[/math] — нормированное векторное пространство данных и [math]\displaystyle{ Y }[/math] — нормированное векторное пространство решений. Для некоторой точки [math]\displaystyle{ x \in X }[/math] задача называется плохо обусловленной, если малое возмущение в [math]\displaystyle{ x }[/math] приводит к большому изменению значения [math]\displaystyle{ f(x) }[/math]. Мы можем количественно определить это, определив число обусловленности, которое показывает, насколько хорошо обусловлена проблема. Определим его как
- [math]\displaystyle{ \widehat{\kappa} = \lim_{\delta \to 0} \sup_{\| \delta x \| \leq \delta} \frac{\| \delta f \|}{\| \delta x \|}. }[/math]
Неустойчивость — это склонность компьютерных алгоритмов, зависящих от арифметики с плавающей запятой, получать результаты, резко отличающееся от точного математического решения задачи. Когда матрица содержит реальные данные со многими значащими цифрами, многие алгоритмы для решения таких задач, как линейные системы уравнений или оптимизация методом наименьших квадратов, могут давать очень неточные результаты. Создание устойчивых алгоритмов для плохо обусловленных задач — центральная задача численной линейной алгебры. Одним из примеров: устойчивость метода отражения Хаусхолдера делает алгоритм надежным методом решения линейных систем, тогда как неустойчивость метода нормальных уравнений для решения задач наименьших квадратов является причиной для выбора методов матричного разложения, таких как использование SVD-разложения. Некоторые методы разложения матриц могут быть неустойчивыми, но имеют простые модификации, делающие их устойчивыми; например, нестабильный метод Грама-Шмидта, который может легко быть изменён для получения устойчивой модификации Грама-Шмидта[1]:140. Другая классическая задача числовой линейной алгебры заключается в том, что метод Гаусса нестабилен, но становится устойчивым с введением опорного элемента.
Итерационные методы
Есть две причины, по которым итерационные алгоритмы — важная часть численной линейной алгебры. Во-первых, многие численные задачи не имеют прямого решения; чтобы найти собственные значения и собственные векторы произвольной матрицы, мы можем использовать только итеративный подход. Во-вторых, безытерационные алгоритмы для произвольной [math]\displaystyle{ m \times m }[/math] матрицы требуют [math]\displaystyle{ O(m^3) }[/math] времени, что неожиданно долго, учитывая, что матрицы содержат только [math]\displaystyle{ m^2 }[/math] чисел. Итерационные подходы могут использовать некоторые особенности матриц, чтобы сократить это время. Например, когда матрица разреженная, итеративный алгоритм может пропустить многие шаги, которые обязательно последуют при прямом подходе, даже если они избыточны.
Ядро многих итерационных методов в численной линейной алгебре — это проекция матрицы на более низкое подпространство Крылова более низкой размерности, что позволяет аппроксимировать характеристики матрицы высокой размерности путем итеративного вычисления эквивалентных характеристик подобных матриц, начиная с пространства низкой размерности и последовательного перехода к более высоким измерениям. Когда [math]\displaystyle{ A }[/math] симметрична, и мы хотим решить линейную задачу [math]\displaystyle{ Ax=b }[/math], классический итеративный подход — метод сопряжённых градиентов. Если [math]\displaystyle{ A }[/math] не симметричная, то примерами итерационных решений линейной задачи являются обобщенный метод минимальных невязок (GMRES)[англ.] и метод сопряжённых градиентов для нормальных уравнений (CGN)[англ.] для нормальных уравнений. Если [math]\displaystyle{ A }[/math] симметричная, то для нахождения собственных значений и собственных векторов мы можем использовать алгоритм Ланцоша[англ.], и если [math]\displaystyle{ A }[/math] не симметричная, то используется итерация Арнольди.
Программное обеспечение для численного анализа
Некоторые языки программирования используют методы оптимизации численной линейной алгебры и предназначены для реализации её алгоритмов. Эти языки включают MATLAB, Analytica, Maple, и Mathematica. Другие языки программирования, которые явно не предназначены для численной линейной алгебры, имеют библиотеки, которые предоставляют процедуры и оптимизацию; в C и Fortran есть пакеты Basic Linear Algebra Subprograms и LAPACK, в python есть библиотека NumPy, и в Perl есть Perl Data Language. Многие команды численной линейной алгебры в R полагаются на фундаментальные библиотеки, такие как LAPACK[5].
Ссылки
- ↑ 1,00 1,01 1,02 1,03 1,04 1,05 1,06 1,07 1,08 1,09 1,10 1,11 1,12 1,13 1,14 1,15 1,16 Trefethen Lloyd, Bau III David. Numerical Linear Algebra (1st edition) (англ.). — Philadelphia: SIAM, 1997. — ISBN 978-0-89871-361-9.
- ↑ 2,0 2,1 2,2 2,3 Golub, Gene H. A History of Modern Numerical Linear Algebra. University of Chicago Statistics Department. Дата обращения: 17 февраля 2019.
- ↑ (1947) «Numerical inverting of matrices of high order». Bulletin of the American Mathematical Society 53 (11): 1021–1099. doi:10.1090/s0002-9904-1947-08909-6.
- ↑ 4,0 4,1 4,2 Golub, Gene H. Matrix Computations / Gene H. Golub, Charles F. Van Loan. — 3rd. — Baltimore : The Johns Hopkins University Press, 1996. — ISBN 0-8018-5413-X.
- ↑ Rickert, Joseph R and Linear Algebra. R-bloggers (August 29, 2013). Дата обращения: 17 февраля 2019.