QR-разложение с возрастающими ошибками
Я получил код для работы, и в нем есть подпрограмма декомпозиции QR, как показано ниже. При проверке этого QR-разложения с использованием расчета ошибок, который я добавил в конце, я обнаружил, что ошибка увеличивается до очень высоких величин после сотен шагов моделирования.
QR-разложение используется для решения коэффициентов, используемых в триквадратичной интерполяции, адаптированной из трилинейной интерполяции, найденной в Википедии: https://en.wikipedia.org/wiki/Trilinear_interpolation. В коде мы генерируем сетку значений 5x5x5 и используем триквадратическую интерполяцию для получения формы пространственной функции значений, чтобы получить f(x,y,z) = a_0 + a_1x + a_2y + a_3*z +...
Функцию f(x,y,z) можно записать как Ax = b для ранее существовавших значений в каждой точке сетки, так что b — вектор существующих значений размером 125x1, A — матрица размером 125x20, рассчитанная с использованием координат, и x — вектор коэффициентов 20x1, состоящий из (a_0, a_1, a_2, ...).
Затем x рассчитывается с использованием алгоритма QR-разложения, показанного ниже. Когда я искал QR-разложение, я нашел здесь некоторые обсуждения, в которых упоминалось, что некоторые подходы могут становиться все более нестабильными, но такие подходы, как размышление домохозяина, должны быть очень стабильными.
Мой первый вопрос: у меня возникли проблемы с определением типа подхода, используемого в моем коде, он не нормализует переменные, и я не могу сопоставить этот процесс с другими реализациями, найденными в Интернете. Может ли кто-нибудь сказать мне, является ли эта реализация уже размышлением домохозяина?
Мой второй вопрос: если реализация правильна и должна быть стабильной, может ли кто-нибудь дать мне несколько идей, чтобы проверить, почему решение будет становиться все более нестабильным? Или, если это не «домохозяин» или просто плохая реализация, поможет ли пример реализации «домохозяин», который я нашел в Интернете, повысить точность?
! QR decomposition solving for vector B in AX = B
subroutine qr_eqsystem(A, B, X, M, N, error)
implicit none
integer, intent(in) :: m, n
real, intent(in) :: A(m, n), B(m)
real, intent(out) :: x(n)
real, intent(out) :: error
!locals
real :: q(m, n), r(n, n), aux
integer :: i, j, k
do i = 1,n
do k = 1,m
q(k, i) = A(k, i)
enddo
enddo
do i = 1,n
do j = 1,n
r(j, i) = 0.0
enddo
enddo
do i = 1,n
do j = 1,i-1
r(j, i) = 0.0
do k = 1,m
r(j, i) = r(j, i) + q(k, j)*A(k, i)
enddo
do k = 1,m
q(k, i) = q(k, i) - r(j, i)*q(k,j)
enddo
enddo
r(i, i) = 0.0
do k = 1,m
r(i, i) = r(i, i) + q(k, i)*q(k, i)
enddo
r(i, i) = sqrt(r(i, i))
if (r(i,i).ne.0.0)then
do k = 1,m
q(k, i) = q(k, i) / r(i, i)
enddo
endif
enddo
do i = 1,n
x(i) = 0.0
do j = 1,m
x(i) = x(i) + q(j, i)*B(j)
enddo
enddo
if (r(n, n).ne.0.0)then
x(n) = x(n) / r(n,n)
endif
do i = n-1,1,-1
do j = i+1,n
x(i) = x(i) - r(i, j) * x(j)
enddo
if (r(i, i).ne.0.0)then
x(i) = x(i) / r(i, i)
endif
enddo
error = 0
! calculate AX and compare with B as an error, report total raw error
do i = 1,m
aux = 0.0
do j = 1,n
aux = aux + A(i,j)*X(j)
enddo
error = error + abs(aux - B(i))
enddo
endsubroutine qr_eqsystem
@MartinBrown, спасибо, Мартин, Real компилируется с двойной точностью, поэтому это не должно быть проблемой, связанной с точностью. Согласно комментарию по подбору полиномов Чебышева, вы имеете в виду подход триквадратичной интерполяции, который использует полиномы типа x^n, и предлагаете использовать синусы и косинусы?
Измените масштаб каждого диапазона переменных на -1,1 для x,y,z, а затем используйте cos(n.acos(x)) вместо x^n. Полученная матрица подгонки будет гораздо лучше обусловлена. На самом деле я никогда не делал этого для нескольких переменных, но это очень хорошо работает для сложных в других отношениях функций в 1D. Впоследствии вы можете преобразовать обратно в полиномиальную форму, если это то, что вам нужно.
«real компилируется с двойной точностью». Если вы используете флаги компилятора для изменения типа переменных, вы запрашиваете огромную боль. Пожалуйста, покажите нам, как именно вы скомпилировали код. И почему вы не используете LAPACK?
Я еще не совсем проснулся, но приведенное выше похоже на подход Грама-Шмидта, который определенно является одним из численно нестабильных подходов. Я настоятельно рекомендую вам использовать LAPACK. Никто не должен был писать плотные процедуры линейной алгебры за последнюю четверть века.
Чтобы сказать гораздо больше, нам нужен минимальный, воспроизводимый пример
Возможно, вы уже используете что-то подобное, но если вы еще этого не видели, это довольно хороший трюк для улучшения состояния матрицы, даже если вы придерживаетесь x^n. Выбор правильных значений x может помочь при интерполяции полиномов.
@IanBush, спасибо, Ян, за идеи. Я рассмотрю возможность использования LAPACK, но с кодом на Фортране, который я использую, это может быть очень сложно. Первоначально код был разработан в 1993 году, и с тех пор над ним работали многие исследователи. Он просто компилируется с флагами «real-size:64» и «Qsave». Мой наиболее вероятный подход — найти реализацию отражения домохозяина и вместо этого добавить ее в код.
@IanBush Я прочитал страницу установки LAPACK (я искал «LAPACK для Windows») и у меня есть несколько вопросов по установке LAPACK в Windows. Не могли бы вы мне с этим помочь? Может быть, нам стоит перейти в личные сообщения?
Извини! Я вообще не использую Windows
@IanBush, нет проблем, ты уже очень помог. Тогда только один вопрос: нужны ли мне все BLAS, LAPACK и LAPACKE?
Если вы, ребята, хотите оставить ответ ниже, я приму его.
Ответы 1
Ваша реализация QR-декомпозиции выглядит нормально (далеко от идеала, но нормально). Вы можете сравнить абсолютные значения диагонали R с реализацией LAPACK или с Python numpy.linalg.qr().
Однако численная устойчивость решения означает лишь то, что
|x - x'|/|x| <= f(m, n)*cond(A)*eps
где f(m, n) — функция малой степени от m и n, cond(A) — число обусловленности, примерно равное отношению диагонали максимального и минимального абсолютных значений матрицы R (A ≈ QR).
Однако величина невязки решения системы наименьших квадратов Ax = B все еще может быть значительной по двум причинам:
- Система несовместима;
- Количество условностей (
cond(A)) слишком велико.
Думаю, стоит поискать вычислительные ошибки при вычислении матрицы A и правильных частей B.
Я думаю, вы правы в том, что проблема, по крайней мере частично, связана с тем, как была вычислена матрица A, и уточнения этой стороны проблемы могут принести дивиденды. Подбор полиномов напрямую приводит к очень плохо обусловленной матричной задаче. ОП также может захотеть попробовать компилятор, который обеспечивает собственную точность x87 или даже четырехкратную точность для его решателя матриц - это значительно уменьшит числовую ошибку в решателе при незначительном изменении кода real -> real*10 или real*16. ISTR Silverfrost (ранее Salford) и GNU имеют такую возможность. Код пытается найти равный полином пульсации?
@MartinBrown, привет, Мартин, я начал новый вопрос, поскольку моя проблема изменилась: stackoverflow.com/questions/78721961/… Я проверил разложенные матрицы Q и R, и оказалось, что само разложение не является проблемой, но каким-то образом обратная замена, решенная x в Ax=B, приводит к большой ошибке при оценке Ax - B
Привет, Серж, пожалуйста, смотрите комментарий выше, чтобы найти ссылку на новый вопрос с дополнительной информацией. Обе предложенные вами причины, несовместимость системы и обусловленность, влияют на сам процесс декомпозиции QR, это правильно? В моем новом вопросе описывалось, что я проверил, что Q и R являются разумными.
Если
realв выбранном вами языке на самом деле не является 64-битным FP с двойной точностью, первое, что нужно попробовать, это увеличить точность массива и рабочих переменных как минимум доdouble, и если ваш язык позволяет это x87, 80-битные числа с плавающей запятой. Кстати, вам, вероятно, лучше подогнать полиномы Чебышева, а не x^n, y^n, z^n. Я думаю, что ваши проблемы могут быть связаны с некорректной постановкой вопроса, который задается решателю матричных уравнений, а не самому решателю. Самый простой способ проверить — найти другой решатель, совместимый по выводам, и попробовать его.