Информационная система "Конференции"

Международная конференция молодых ученых по математическому моделированию и информационным технологиям

29-31 октября 2002 г., Новосибирск, Академгородок

Мацех А.М.
ИВТ СО РАН (Новосибирск)

Анализ алгоритмов вычисления сингулярных чисел и сингулярных векторов несимметричных матриц по соответствующим матрицам в форме Голуба-Кахана.

Аннотация:

Сингулярные числа произвольной несимметричной матрицы являются еe ортогональными инвариантами и совпадают с сингулярными числами некоторой двухдиагональной матрицы, ортогонально подобной исходной. Сингулярные числа произвольной двухдиагональной матрицы в свою очередь совпадают с неотрицательными собственными значениями симметричной трeхдиагональной матрицы в так называемой форме Голуба-Кахана, имеющей нулевую главную диагональ и побочные диагонали составленные из перемежевающихся элементов соответствующей двухдиагональной матрицы. В работе исследуются вопросы адаптации алгоритмов решения симметричных задач на собственные значения с гарантированной точностью для определения сингулярных чисел и сингулярных векторов несимметричных матриц по соответствующим матрицам в форме Голуба-Кахана.

Сингулярное разложение прямоугольных матриц

Произвольная прямоугольная вещественная матрица $ A \in R^{N \times M}$ может быть преведена к спектрально-эквивалентному двухдиагональному виду

$\displaystyle \tilde{D} = P^{*} A\, Q $

некоторыми ортогональными преобразованиями $ P$ и $ Q$ [Год97]. Ненулевая часть $ D$ матрицы $ \tilde{D}$ является квадратной двухдиагональной матрицей размерности $ n \times n, \ n = min\{N, M\}$. Без потери общности можно предположить что матрица $ D$ является верхнедиагональной матрицей вида

$\displaystyle D = \left[\begin{array}{ccccc} d_1 & b_1 & & & \\ & d_2 & b_2 &... ... & \ddots & \\ & & & d_{n-1} & b_{n-1} \\ & & & & d_n \end{array} \right], $

причем элементы главной диагонали $ d_1, d_2, \ldots, d_n$ являются собственными значениями матрицы $ D$. Произвольная квадратная вещественная матрица допускает сингулярное разложение [Год97], и следовательно матрица $ D$ может быть представлена в виде

$\displaystyle D = V\, \Sigma\, U^{*}, $

где $ \Sigma$ является диагональной матрицей вида

$\displaystyle \Sigma = \left[\begin{array}{ccccc} \sigma_1 & & & \\ & \sigma_2 & & \\ & &\ddots& & \\ & & & & \sigma_{n} \\ \end{array} \right] $

составленной из сингулярных чисел $ \sigma_i \geq 0, \ i=1, 2, \ldots, n,$ а ортогональные матрицы $ U$ и $ V$ называются соответственно матрицами левых и правых сингулярных векторов. Таким образом сингулярное разложение исходной матрицы $ A$ имеет вид

$\displaystyle A = (P\,V)\, \Sigma\, (Q\,U)^{*}, $

что означает что задача об определении сингулярного разложения для матрицы $ A$ сводится к задаче определения сингулярного разложения для двухдиагональной матрицы $ D$ [ГАКК88], которая может быть получена к примеру применением цепочек ортогональных преобразований отражения если матрица $ A$ плотная, и методом Арнольди или несимметричным методом Ланцоша, если матрица $ A$ является разреженной матрицей с известной структурой [GL96]. В этой работе мы не останавливаемся на деталях реализации этих алгоритмов, концентрируя внимание на решении соответствующей двухдиагональной задачи. Сингулярные числа этих задач совпадают с точностью определяемой качеством ортогональных преобразований $ P$ и $ Q$. В общем погрешность таких преобразованой должна быть величиной порядка погрешности машинного представления единицы умноженной на норму матрицы $ A$. Применяя преобразования $ P$ и $ Q$ соответственно к левым и правым сингулярным векторам матрицы $ D$ мы получаем аппроксимацию сунгулярных векторов матрицы $ A$, погрешность которых уже зависит не только от погрешностей двухдиагонализации, но также и от точности определения сингулярного разложения матрицы $ D$.

Связь между сингулярным и собственными числами

Сингулярные числа произвольной матрицы могут быть определены как неотрицательные собственные значения соответствующей формы Жордана-Вейланда [Ste93]. Форма Жордана-Вейланда для двухдиагональной матрицы $ D$ имеет вид

$\displaystyle J = \left[ \begin{array}{cc} 0 & D \\ D^* & 0 \end{array} \right] $

причем собственные значения $ \lambda(J)$ и сингулярные числа $ \sigma(D)$ связяны следующим образом:

$\displaystyle \lambda(J) = \{ \sigma(D) \} \cup \{ -\sigma(D) \}. $

Ортогональными преобразованиями перестановок $ \Pi$ форма Жордана-Вейланда $ J$ двухдиагональной матрицы $ D$ может быть приведена к эквивалентной форме Голуба-Кахана [Fer98]

$\displaystyle G = \left[\begin{array}{cccccc} 0 & d_1 & & & & \\ d_1 & 0 & b_... ...ots & & \\ & & & b_{n-1} & 0 & d_n \\ & & & & d_n & 0 \end{array} \right] $

причем связь между собственными значениями и сингулярными числами сохраняется:

$\displaystyle \lambda(G) = \{ \sigma(D) \} \cup \{ -\sigma(D) \}, $

a собственный вектор $ (v, u)^{T}$ матрицы $ J$ переводится преобразованием $ \Pi$ в вектор [ГАКК88]

$\displaystyle \left[\begin{array}{c} u^1 \\ v^1 \\ u^2 \\ v^2\\ \vdots \\ ... ...\end{array} \right] = \Pi \left[\begin{array}{c} v \\ u \end{array} \right]. $

Следовательно, решая частичную задачу на собственные значения для матрицы $ G$ мы можем определить сингулярные числа и соответствующие сингулярные векторы матрицы $ D$.

Вычисление сингулярного разложения матрицы Голуба-Кахана

Вычислительная схема сингулярного разложения произвольной вещественной матрицы $ A$ по соответствующей матрице Голуба-Кахана может быть сформулирована в виде следующего алгоритма:

Алгоритм 1 (Сингулярное разложение)  
  1. Определить ортогональные преобразования $ P$ и $ Q$ позволяющие представить прямоугольную матрицу $ A \in R^{N \times M}$ в эквивалентной двухдиагональной форме $ \tilde{D} = P^{*}\, A\, Q $ с ненулевой частью $ D, {\rm dim}(D) = n \times n, \ n = min\{N, M\}$.
  2. Определить $ n$ наибольших собственных значений $ \lambda_{n}, \lambda_{n+1}, \ldots, \lambda_{2 n}$ и соответствующих им собственных векторов $ \omega_{n}, \omega_{n+1}, \ldots, \omega_{2 n}$ трехдиагональной матрицы Голуба-Кахана $ G, {\rm dim}(G) = 2n \times 2 n$, соответствующей матрице $ D$.
  3. Взять собственные значения $ \lambda_{n}, \lambda_{n+1}, \ldots, \lambda_{2 n}$ матрицы $ G$ в качестве сингулярных чисел матрицы $ A$: $ \sigma_1(A) = \lambda_{n}(G),\sigma_2(A) = \lambda_{n+1}(G), \ldots, \sigma_n(A) = \lambda_{2 n}(G)$.
  4. Определить левые сингулярные векторы $ u_k, k = 1, \ldots, n$ матрицы $ D$ через нечетные компоненты собственных векторов $ w_k, k = 1, \ldots, n$ матрицы $ G$, и правые сингулярные векторы $ v_k, k = 1, \ldots, n$ матрицы $ D$ через четные компоненты собственных векторов $ w_k, k = 1, \ldots, n$ матрицы $ G$.
  5. Определить левые сингулярные векторы $ S$ матрицы $ A$ действуя преобразованием $ P$ слева на сингулярные векторы $ U$ матрицы $ D$: $ S = P\, U$ и правые сингулярные векторы $ Т$ матрицы $ A$ действуя преобразованием $ Q$ слева на сингулярные векторы $ V$ матрицы $ D$: $ S = Q\, V$.

Остановимся подробнее на деталях машинной реализации описанной схемы. Мы предполагаем что матрица $ A$ приводится к двухдиагональному виду $ D$ почти ортогональными преобразованиями с погрешностью не намного превышающей $ \epsilon \Vert A\Vert$, где $ \epsilon$ представляет собой погрешность округления чисел с вещественной точкой на используемой архитектуре. Матрица $ D$ в явном виде не строятся - элементы $ D$ хранятся в виде вектора $ g = (d_1, b_1, d_2, \ldots, b_{n-1}, d_n)$ составленного из перемежевающихся элементов главной и побочной диагоналей матрицы $ D$. Вектор $ g$ соответствует побочной диагонали матрицы Голуба-Кахана $ G$. Сингулярные числа матрицы $ D$ вычисляются методом бисекций основанном на применении закона инерции Сильвестра для подсчета неположительных собственных значений матрицы $ G - I \mu$. Точность при подсчете собственных значений имеет критическое значение при решении данной задачи, поскольку даже малые погрешности при вычислении сингулярных (собственных) чисел влекут за собой значительные погрешности при определении сингулярных (собственных) векторов. Метод бисекций является наиболее робастным среди методов определения собственных чисел симметричных матриц, и при аккуратном подсчете смен знаков позволяет вычислять интервалы содержащие собственные числа рассматриваемой матрицы с гарантированой точностью. Соответствующие собственные векторы матрицы $ G$ могут быть вычислены методом Обратной Итерации [ABB+95], методом Годунова, основанном на применении составных последовательностей Штурма [ГАКК88], или методом Годунова с итерационным уточнением [МШ01].

Особенности вычисления собственных значений матрицы Голуба-Кахана

Традиционно подсчет неотрицательных собственных чисел основывался на подсчете неположительных элементов в последовательности Штурма матрицы $ G - I \mu$ [Год97]. Однако, поскольку последовательности Штурма вычисляются по главным минорам матрицы $ G - I \mu$, среди элементов последовательности могут появляться числа непредставимые в машинной арифметике. Для избежания таких ситуаций при подсчете собственных чисел по последовательностям Штурма требуется вносить малые погрешности в элементы матрицы $ G$ [ГАКК88], имеющей нулевую главную диагональ. Однако внесение погрешностей в элементы главной диагонали нежелательно. Современная версия метода бисекций основывается на применении закона инерции Сильвестра, согласно которому инерция конгруэнтных матриц остается неизменной [GL96]. Это означает что число отрицательных, положительных и нулевых собственных значений симметричной матрицы $ G - I \mu \in R^{n \times n}$ и конгруэнтной ей матрицы $ L^{T}(G - I\mu)L$, где матрица $ L \in R^{n \times n}$ невырожденная, совпадает. Диагональная матрица $ L^{T}(G - I\mu)L$ может быть получена при помощи Гауссова исключения, причем элементы $ q_i$ главной диагонали матрицы $ L^{T}(G - I\mu)L$ имеют вид [Fer98]:

$\displaystyle q_i = - \mu - g_{i-1}^2/q_{i-1}, \ i=2, \ldots, 2n, \ q_1 = -\mu $

В IEEE-арифметике, при правильном масштабировании матрицы $ D$ для избежания ситуаций деления на ноль и бесконечность, функция $ \tau(\mu)$, равная числу отрицательных элементов $ q_i$, является монотонной функцией, и следовательно пригодна для подсчета отрицательных собственных значений матрицы $ (G - I\mu)$  [Fer98], причем элементы матрицы $ D$ масштабируются так, чтобы удовлетворялось условие

$\displaystyle 2^m \max(\max_i d_i, \max_i b_i) \leq \{\gamma/\eta\}^{1/2} $

где $ \eta$ представляет собой минимальное положительное машинное число представимое в IEEE-арифметике, $ m$ является максимальным целым машинным числом, а $ \gamma$ определяется таким образом, что $ \epsilon > \gamma = 2^l \geq \eta$ для некоторого целого числа $ l$. Посколку $ m$ является целым числом, параметр масштабирования $ 2^m$ не вносит погрешностей в элементы матрицы $ D$. Требуя дополнительно, чтобы $ q_i = -\gamma$ в случае если $ -\gamma \leq q_i < \gamma$, мы получаем алгоритм для подсчета отрицательных элементов $ q_i$ матрицы $ L^{T}(G - I\mu)L$:

Алгоритм 2 (Инерция матрицы Голуба-Кахана)  
$ $
$ \tau =$   inertia$ (\mu,\ g)$
$ $
масштабировать $ D$ так, что $ 2^m \max g_i \leq \{\gamma/\eta\}^{1/2}$
$ q = - \mu, \ \tau = 1$
for
$ i = 1, 2, \ldots, 2n-1$
if
$ -\gamma \leq q < \gamma$
then
$ q = -\gamma$
$ q = - \mu - g_{i}^2/q$
if
$ q < 0 $
then
$ \tau = \tau + 1$
end

Сингулярные числа матрицы $ D$ ищутся на интерале $ [X, Y]$, определяемом свойствами сингулярных чисел трехдиагональных матриц следующим образом [ГАКК88]:

$\displaystyle Y = \mathfrak{M}(D),\ X = \vert d_1 d_2\ldots d_n\vert/Y^{n-1}, $

причем

$\displaystyle \mathfrak{M}(D) = \max(\underset{1 \leq i < n}{\max}(d_i + b_i), \underset{1\leq i < n}{\max}(d_{i+1} + b_i)) $

Положив $ \mu^{\prime} = (X+Y)/2$, и вычислив инерцию матрицы $ G$, можно выяснить число собственных значений, меньших $ \mu$, а значит решить какому из полуинтервалов $ [X, \mu^{\prime)}$, или $ (\mu^{\prime}, Y]$, принадлежит $ \mu$. Процесс последовательного деления интервалов, окружающих $ \mu$, продолжается до пор, пока в результате сближения границ величина интервала $ (x_{i}, y_{i})$, содержащего $ \mu_{i}$, не сравнима с $ \epsilon(\vert x_{i}\vert+\vert y_{i}\vert)$.

Вычисление сингулярных векторов.

Поскольку задача отыскания сингулярных векторов двухдиагональной матрицы $ D$ сводится к задаче определения собственных векторов симметричной трехдиагональной матрицы $ G$, мы исследовали возможность преименения методов Обратной Итерации [ABB+95], метода Годунова [ГАКК88], и метода Годунова с итерационным уточнением [МШ01] для решения данной задачи. Метод Годунова имеет сложность $ O(n)$ операций с вещественной точкой для определения собственного вектора, но требует внесения погрешностей в элемены главной диагонали матрицы $ D$, что приводит к потере точности при вычислениях. Одновременно метод Годунова генерирует почти коллинеарные векторы, соответствующие близким собственным значениям, в то время как метод Обратной Итерации включает в себя шаг выборочной реортогонализации. В большенстве экспериментов метод обратной итерации доставляет меньшие погрешности, чем метод Годунова, являющийся прямым методом по своей природе, и требующим внесения погрешностей в нулевые диагональные элементы. Однако метод Обратной Итерации является итерационной процедурой требующей вплоть до $ O(n^2)$ действий с вещественной точкой для вычисления собственного вектора. В отличии от алгоритма Годунова в методе Обратной Итерации сходимость к решению не гарантирована, поскольку начальный вектор может не иметь нетривиальной компоненты в направлении искомого решения.

Альтернативным подходом к решению данной задачи является разработанный нами гибридный алгоритм Годунова с итерационным уточнением, или метод Годунова-Обратной Итерации [ГАКК88]. Идея метода заключается в том, чтобы использовать векторы расчитанные методом Годунова с гарантированной точностью как черезвычайно точные начальные векторы в методе Обратной Итерации, имеющие нетривиальную компоненту в направлении искомого решения. Как правило одного-двух шагов обратной итерации оказывается достаточным для достижения требуемой точности, в то время как, к примеру метод Обратной Итерации, реализованный в пакете LAPACK в виде процедуры STEIN [ABB+95] требует трех-пяти шагов. Одновременно, на шаге итерационного уточнения мы реортогонализуем векторы вычисленные методом Годунова, получая таким образом результаты сравнимые или даже превосходящие результаты обратной итерации, требующей большего числа шагов, и не всегда сходящегося. Формально метод Годунова-Обратной Итерации может быть описан в виде следующего алгоритма:

Алгоритм 3 (Метод Годунова с Итерационным Уточнением)   Вычисление собственных векторов $ u_{i},\ i=1,\ldots,n$ трехдиагональной матрицы $ T = T^{T} \in R^{n \times n}$ имеющей главную диагональ $ d$ и побочную диагональ $ b$ на архитектуре с погрешностью округления $ \epsilon_{mach}$ и $ t$ битами мантиссы в представлении машинных чисел.
godunov_inverse_iteration(d, b, i, n)
$ $
определить собственные числа $ \mu_{i} \in (\alpha_{i}, \beta_{i})$ методом бисекций так что
$ \vert\beta_{i} - \alpha_{i}\vert\leq \epsilon_{mach}\ \Vert T\Vert _{2}, \ i=1, \ldots, n$
$ k=0, \ \delta = 2^t/(100*n)$
for
$ (i=1, i<=n, i++)$
определить левосторонную последовательность Штурма $ P^{+}_{n}(\alpha_{i})$
$ P^{+}_{0}(\alpha_{i}) = \vert b_{0}\vert / (d_{0} - \alpha_{i})$
$ P^{+}_{k}(\alpha_{i}) =\vert b_{k}\vert/(d_{k}-\alpha_{i}-\vert b_{k-1}\vert) P^{+}_{k-1}(\alpha_{i})$
$ P^{+}_{n-1}(\alpha_{i}) = 1/(d_{n-1} - \alpha_{i} - \vert b_{n-2}\vert) P^{+}_{n-2}(\alpha_{i})$
определить правосторонную последовательность Штурма $ P^{-}_{n}(\beta_{i})$
$ P^{-}_{n-1}(\beta_{i}) = d_{n-1} - \beta_{i}$
$ P^{-}_{k}(\beta_{i}) =(d_{k}-\beta_{i}-\vert b_{k}\vert/P^{-}_{k+1}(\beta_{i}))/\vert b_{k-1}\vert$
$ P^{-}_{0}(\beta_{i}) = d_{0} - \beta_{i} - \vert b_{0}\vert/P^{-}_{1}(\beta_{i})$
определить целое число $ l$, такое, что
$ (P_l^{+}(\alpha_i) - P_l^{-}(\beta_i)) (1/P_{l+1}^{-}(\beta_i) - 1/P_{l+1}^{+}(\alpha_i)) \leq 0$
определить двусторонную последовательность Штурма $ P_{n}(\mu_i)$
$ P_{0}(\mu_{i}), \dots, P_{n-1}(\mu_{i}) \stackrel{\rm def}{=} P^{+}_{0}(\alpha... ...ots P^{+}_{l}(\alpha_{i}), P^{-}_{l+1}(\beta_{i}),\ldots P^{-}_{n-1}(\beta_{i})$
вычислить собственный вектор $ u^{i}$: $ u^{i}_{0}=1,\ u^{i}_{k} =- u^{i}_{k-1}$   sign$ (b_{i-1}) / P_{i-1}(\mu_i)$
$ \gamma_{i} = \beta_{i} $
if
$ (k > 0 \cap \vert\gamma_{i} - \gamma_{i-1}\vert \leq 10\ \epsilon_{mach}\ \vert\gamma_{i}\vert)$
then
$ \gamma_{i}=\gamma_{i-1} + 10\epsilon_{mach}\vert\gamma_{i}\vert$
do
$ $
$ k = k+1$
Решить $ (T - \gamma_{i} I) z^{k} = u_{i}^{k-1}$
for
$ (j=1, j <= k, k++)$
if
$ \vert\gamma_j - \gamma_k\vert \leq \Vert T\Vert _{\infty}/1000 $
then
$ z_{k}=z_{k}-(z_{k}, z_{j}) z_{j}$
end
$ u_{i}^{k} = z^{k}/\Vert z^{k}\Vert _{2}$
while
$ ( \Vert z^{k}\Vert _{\infty} > \delta )$
end

Примеры расчета сингулярных чисел и сингулярных векторов двухдиагональных матриц

Нами разработан пакет программ реализующий описанные алгоритмы сингулярного разложения в двойной точности на языке C в IEEE арифметике. В качестве иллюстрации работы этого пакета мы приводим примеры расчета сингулярных чисел и сингулярных векторов нескольких двухдиагональные матриц (примеры 1- 4), описанных в монографии С. К. Годунаова et. al [ГАКК88]. Расчеты проводились в операционной системе Linux, с использованием компилятора GNU C на персоональном компьютере имеющем следующие параметры: i686 Intel(R) Xeon(TM) CPU 1500MHz.

В таблице 1 представлены погрешности расчета машинных интервалов $ [x_i, y_i]$ содержащих сингулярные числа $ \sigma_i $ с гарантированной точностью, а также максимальные и минимальные сингулярные числа рассматриваемых тестовых матриц. Сингулярные числа $ (\sigma_1, \ldots, \sigma_n)$ были расчитаны методом бисекций, основанном на использовании инерции Сильвестра. В таблице 2 представлены оценки вычислительной погрешности сингулярного разложения $ \Vert A_k V - U \Sigma\Vert _{\infty} $ тестовых матриц, а также оценки ортогональности сингулярных векторов $ U, \ V$ тестовых матриц, расчитанных в соответствии с методом Годунова с Итерационным Уточнением с использованием в качестве сдвига правых границ $ y_i$ интервалов $ [x_i, y_i]$, содержащих соответствующие собственные значения. В таблице 3 представлены алалогичные характеристики сингулярных векторов, расчитанных методом обратной итерации по алгоритму STEIN из пакета LAPACK [ABB+95] с использованием в качестве сдвига правых границ $ y_i$ интервалов $ [x_i, y_i]$, содержащих соответствующие собственные значения, и наконец в таблице  4 представлены соответствующие характеристики сингулярных векторов расчитанных методом Годунова. Сравнивая таблицы 2- 4 мы убеждаемся что метод Годунова с Итерационным Уточнением значительно улучшает качество сингулярных (собственных) векторов, расчитанных методом Годунова. Используя меньшее число итераций чем метод Обратной Итерации, метод Годунова с Итерационным Уточнением позволяет расчитывать сингулярные (собственные) векторы практически с такой же точностью.

Пример 1   % latex2html id marker 2857 $ U^T A_{\ref{ex_a}}\, V,\ {\rm dim}(A_{\ref{ex_a}}) = n \times n,\ n = 1000.$

% latex2html id marker 2859 $\displaystyle A_{\ref{ex_a}} = \left[\begin{array}{... ...\ & &\ddots & \ddots & \\ & & & 1 & 10 \\ & & & & 1 \end{array} \right], $

Сингуярные числа рассматриваемой матрицы были расчитаны с гарантированной точностью $ 7.9936057773011271e-15 $. Результаты аналогичного расчета сингулярных чисел в пакете Matlab согласовались с нашими результатами с максимальной погрешностью равной $ 1.4210854715202004e-14$. Максимальная погрешность сингулярного разложения вычисленного методом Годунова с Итерационным Уточнением равна $ 1.4665441894424352e-14$.

Пример 2   % latex2html id marker 2868 $ U^T A_{\ref{ex_b}}\, V,\ {\rm dim}(A_{\ref{ex_b}}) = n \times n,\ n = 1000.$

% latex2html id marker 2870 $\displaystyle A_{\ref{ex_b}} = \left[\begin{array}{... ...\ddots & \ddots & \\ & & & 0.01 & 900 \\ & & & & 0.01 \end{array} \right], $

Сингуярные числа рассматриваемой матрицы были расчитаны с гарантированной точностью $ 6.2527760746888816e-13$. Результаты аналогичного расчета сингулярных чисел в пакете Matlab согласовались с нашими результатами с максимальной погрешностью равной $ 5.2295945351943374e-12$. Максимальная погрешность сингулярного разложения вычисленного методом Годунова с Итерационным Уточнением равна $ 1.8878383748519391e-12$.

Пример 3   % latex2html id marker 2879 $ U^T A_{\ref{ex_c}}\, V,\ {\rm dim}(A_{\ref{ex_c}}) = n \times n,\ n = 1000.$

% latex2html id marker 2881 $\displaystyle A_{\ref{ex_c}} = \left[\begin{array}{... ... &\ddots & \ddots & \\ & & & 0.5 & 0.5 \\ & & & & 0.5 \end{array} \right], $

Сингулярные числа матрицы % latex2html id marker 2883 $ A_{\ref{ex_c}}$ являются корнями полиномов Чебышева 2-го рода, и выражаются формулой:

$\displaystyle \sigma_i = \frac{\cos(n-i)\, \pi}{2n+1}, \ i = 0, 1, \ldots, n-1. $

Сингуярные числа рассматриваемой матрицы были расчитаны с гарантированной точностью $ 6.9388939039072284e-16$, причем максимальное отклонение вычисленных нами сингулярных чисел от соответствующих значений расчитанных по аналитической формуле равнялось точно $ 2.2204460492503131e-16$ -- т.е. машинной погрешности представления единицы в IEEE арифметике с двойной точностью. Результаты аналогичного расчета сингулярных чисел в пакете Matlab согласовались с аналитическими результатами с максимальной погрешностью равной лишь $ 6.6613381477509392e-16$. Максимальная погрешность сингулярного разложения вычисленного методом Годунова с Итерационным Уточнением равна $ 1.6391144039830456e-15$.

Пример 4   % latex2html id marker 2896 $ U^T A_{\ref{ex_d}}\, V,\ {\rm dim}(A_{\ref{ex_d}}) = n \times n,\ n = 1000.$

% latex2html id marker 2898 $\displaystyle A_{\ref{ex_d}} = \left[\begin{array}{... ... & \ddots & \\ & & & d_{n-1} & b_{n-1} \\ & & & & d_n \end{array} \right], $

Элементы гавной и побочной диагоналей этой матрицы вычисляются по формулам:

$\displaystyle d_i = (i+1)/\sqrt{(2i+1)(2i+3)},\ i = 0, 2, \ldots, 2n-2, $

$\displaystyle b_i = (i+1)/\sqrt{(2i+1)(2i+3)},\ i = 1, 3, \ldots, 2n-3, $

Сингулярные числа матрицы % latex2html id marker 2904 $ A_{\ref{ex_d}}$ являются ортонормированными полиномами Лежандра и были расчитаны с гарантированной точностью $ 7.9936057773011271e-15 $. Результаты аналогичного расчета сингулярных чисел в пакете Matlab согласовались с нашими результатами с максимальной погрешностью равной $ 1.1102230246251565e-15 $. Максимальная погрешность сингулярного разложения вычисленного методом Годунова с Итерационным Уточнением равна $ 2.1312193757624329e-15$.


Таблица: Погрешность определения сингулярных чисел, и максимальные и минимальные сингулярные числа матриц из примеров % latex2html id marker 2972 $ {\ref{ex_a}}-{\ref{ex_d}}$.
матрица $ \underset{1\leq i \leq n}{\max}([x_i, y_i] \ni \sigma_i)$ $ \sigma_{max} $ $ \sigma_{min} $
% latex2html id marker 2918 $ A_{\ref{ex_a}} $ $ 7.9936057773011271e-15 $ $ 1.0999995514634513e+01 $ $ 1.4300021896728732e-18 $
% latex2html id marker 2926 $ A_{\ref{ex_b}} $ $ 6.2527760746888816e-13$ $ 9.0000999995065263e+02 $ $ 3.3069079079233366e-15 $
% latex2html id marker 2934 $ A_{\ref{ex_c}} $ $ 6.9388939039072284e-16$ $ 9.9999876753247885e-01 $ $ 7.8500557994265236e-04 $
% latex2html id marker 2942 $ A_{\ref{ex_d}} $ $ 7.2164496600635175e-16 $ $ 9.9999927746317030e-01 $ $ 7.8520175772144778e-04 $



Таблица: Вычислительная погрешность сингулярного разложения и характеристика ортогональности левых и правых сингулярных векторов матриц из примеров % latex2html id marker 3037 $ {\ref{ex_a}}-{\ref{ex_d}}$ вычисленных методом Годунова с Итерационным Уточнением.
матрица $ \Vert A_k V - U \Sigma\Vert _{\infty} $ $ \Vert V^T V - I\Vert _{\infty} $ $ \Vert U^T U - I\Vert _{\infty} $
% latex2html id marker 2983 $ A_{\ref{ex_a}} $ $ 1.4665441894424352e-14$ $ 1.0421333544105856e-14 $ $ 1.0360716101911912e-14 $
% latex2html id marker 2991 $ A_{\ref{ex_b}} $ $ 1.8878383748519391e-12$ $ 2.7646697944860962e-15 $ $ 2.9595438227334469e-15 $
% latex2html id marker 2999 $ A_{\ref{ex_c}} $ $ 1.6391144039830456e-15$ $ 6.4705650678520278e-15 $ $ 6.5167603348807169e-15 $
% latex2html id marker 3007 $ A_{\ref{ex_d}} $ $ 2.1312193757624329e-15$ $ 5.5689836596815914e-15 $ $ 5.4758065424539648e-15 $



Таблица: Вычислительная погрешность сингулярного разложения и характеристика ортогональности левых и правых сингулярных векторов матриц из примеров % latex2html id marker 3102 $ {\ref{ex_a}}-{\ref{ex_d}}$ вычисленных методом Обратной Итерации.
матрица $ \Vert A_k V - U \Sigma\Vert _{\infty} $ $ \Vert V^T V - I\Vert _{\infty} $ $ \Vert U^T U - I\Vert _{\infty} $
% latex2html id marker 3048 $ A_{\ref{ex_a}} $ $ 1.9194750185848903e-14 $ $ 9.6737639915873590e-15 $ $ 9.6762785615362646e-15 $
% latex2html id marker 3056 $ A_{\ref{ex_b}} $ $ 1.5694406873168894e-12 $ $ 2.5867578327518879e-15 $ $ 2.8761806369179954e-15 $
% latex2html id marker 3064 $ A_{\ref{ex_c}} $ $ 1.4384580160596148e-15 $ $ 6.4933559252047057e-15 $ $ 6.5132519654545167e-15 $
% latex2html id marker 3072 $ A_{\ref{ex_d}} $ $ 1.5327610597247947e-15 $ $ 5.8238706503085389e-15 $ $ 5.7114242318572437e-15 $



Таблица: Вычислительная погрешность сингулярного разложения и характеристика ортогональности левых и правых сингулярных векторов матриц из примеров % latex2html id marker 3167 $ {\ref{ex_a}}-{\ref{ex_d}}$ вычисленных методом Годунова.
матрица $ \Vert A_k V - U \Sigma\Vert _{\infty} $ $ \Vert V^T V - I\Vert _{\infty} $ $ \Vert U^T U - I\Vert _{\infty} $
% latex2html id marker 3113 $ A_{\ref{ex_a}} $ $ 9.9498743710661994e+00 $ $ 1.0000000000001692e+00 $ $ 4.7912138738717497e-10 $
% latex2html id marker 3121 $ A_{\ref{ex_b}} $ $ 8.9999999994444431e+02 $ $ 1.0000000001277405e+00 $ $ 3.8780738431827434e-06 $
% latex2html id marker 3129 $ A_{\ref{ex_c}} $ $ 8.2748107998609149e-12 $ $ 1.9698795341085592e-10 $ $ 1.9733465285846337e-10 $
% latex2html id marker 3137 $ A_{\ref{ex_d}} $ $ 3.4417467212941555e-13 $ $ 8.5444802200893734e-11 $ $ 8.5449407679558367e-11 $


Список Литературы

ГАКК88
С. К. Годунов, А. Г. Антонов, О. П. Кирилюк, В. И. Костин.
Гарантированная точность решения систем линейных уравнений в евклидовых пространствах.
Наука, Новосибирск, 1988.

Год97
С. К. Годунов.
Современные аспекты линейной алгебры.
Научная книга, Новосибирск, 1997.

МШ01
А. М. Мацех, Э. П. Шурина.
Решение трехдиагональной симметричной проблемы собственных значений с гарантированной точностью за О(n**2) операций.
Современные проблемы прикладной математики и механики: теория, эксперимент и практика. Академгородок, Новосибирск, Россия,  24-29 Июня 2001.

ABB+95
E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A McKenney, S. Ostrouchov, D. Sorensen.
LAPACK Users' Guide.
SIAM, Philadelphia, second edition, 1995.

Fer98
K. V. Fernando.
Accurately counting singular values of bidiagonal matrices and eigenvalues of skew-symmetric tridiagonal matrices.
SIAM Journal on Matrix Analysis and Applications, 20(2):373-399, 1998.

GL96
Gene H. Golub, Charles F. Van Loan.
Matrix Computations.
The Johns Hopkins University Press, third edition, 1996.

Ste93
G. W. Stewart.
On the early history of the singular value decomposition.
SIAM Review, 35(4):551-566, 1993.

About this document ...

Анализ алгоритмов вычисления сингулярных чисел и сингулярных векторов несимметричных матриц по соответствующим матрицам в форме Голуба-Кахана.

This document was generated using the LaTeX2HTML translator Version 2K.1beta (1.47)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -no_math -html_version 3.2 -split 0 ysci_2002.tex ysci_2002.html

The translation was initiated by Anna Matsekh on 2002-10-23

Anna Matsekh 2002-10-23