Почему плохо обусловленные линейные системы могут быть решены точно?

13

Согласно ответу здесь , большое число условий (для линейного решения системы) уменьшает гарантированное количество правильных цифр в решении с плавающей запятой. Матрицы дифференцирования более высокого порядка в псевдоспектральных методах обычно очень плохо обусловлены. Почему же они все еще очень точные методы?

Я понимаю, что низкая точность, исходящая из плохо обусловленных матриц, является только гарантированным значением, но все же меня удивляет, почему плохо обусловленные матрицы точно решаются прямыми методами на практике - например, LCOLстолбцы таблицы 3.1 на стр. 11 Wang et al., ХОРОШИЙ МЕТОД КОЛЛОКАЦИИ С ИСПОЛЬЗОВАНИЕМ МАТРИЦЫ ПСЕВДОСПЕКТРАЛЬНОЙ ИНТЕГРАЦИИ , SIAM J. Sci. Comput., 36 (3) .

Золтан Цати
источник
2
Моя интуиция заключается в том, что разрешимость / точность системы Ax = b связана с вектором воздействия b, а не только с матрицей A. Возможно, если b не "прощупывает" или "возбуждает" плохо обусловленные моды A, то точное решение остается возможным В качестве ограничивающего примера, A может быть в точности единичным (бесконечное число условий), однако Ax = b все еще может иметь решение, которое можно точно вычислить, если форсирующие данные b находятся в диапазоне A. Я допускаю, что это довольно неплохо волнистые, поэтому я только комментирую вместо ответа.
rchilton1980
@ rchilton1980 "но Ax = b может все еще иметь решение" Но это решение не уникально. И примеры, на которые я ссылаюсь, обладают уникальным решением.
Золтан Ксати
Это справедливый контрапункт - возможно, артефакт выбора бесконечного числа условий (ровно нулевого собственного значения). Однако я думаю, что вы можете заменить это нулевое собственное значение машинным эпсилоном, и моя точка зрения остается в силе. (То есть система имеет очень большое число условий, система неособа с уникальным решением, которое мы можем очень точно вычислить, если b не имеет компонента вдоль этой крошечной собственной пары).
rchilton1980
1
Чтобы быть более конкретным, мой мысленный эксперимент здесь выглядит примерно так: A = diag ([1 1 1 1 1 eps]), b = [b1 b2 b3 b4 b5 0]. Он придуман, но я думаю, что достаточно обосновать первоначальное утверждение: «иногда плохо обусловленные А могут быть точно решены для конкретного выбора b»
rchilton1980
1
Просто приведите еще один пример из блога Мёлера blogs.mathworks.com/cleve/2015/02/16/…
percusse

Ответы:

7

Добавлено после моего первоначального ответа:

Теперь мне кажется, что автор упомянутой статьи дает номера условий (очевидно, 2-нормальные номера условий, но, возможно, номера условий бесконечности норм) в таблице, в то же время давая максимальные абсолютные ошибки, а не относительные ошибки нормы или максимальные относительные ошибки по элементам все это разные меры.) Обратите внимание, что максимальная поэлементная относительная ошибка - это не то же самое, что относительная ошибка в бесконечности. Кроме того, ошибки в таблице связаны с точным решением краевой задачи исходного дифференциального уравнения, а не с дискретной линейной системой уравнений. Таким образом, информация, представленная в документе, действительно не подходит для использования с ошибкой, связанной с номером условия.

Однако в моей репликации вычислений я вижу ситуации, когда относительная ошибка нормы бесконечности (или относительная ошибка двух норм) существенно меньше границы, установленной номером условия бесконечности нормы (соответственно, 2-нормальным числом условия). Иногда тебе просто везет.

Я использовал пакет DMSUITE MATLAB и решил примерную задачу из этой статьи, используя псевдоспектральный метод с полиномами Чебышева. Мои показатели состояния и максимальные абсолютные ошибки были аналогичны тем, о которых сообщалось в статье.

Я также видел нормальные относительные ошибки, которые были несколько лучше, чем можно было бы ожидать на основе номера условия. Например, на примере задачи с , используя N = 1024 , я получаюϵ=0.01N=1024

Cond (А, 2) = 7.9e + 8

Cond (А, инф) = 7.8e + 8

норма (U-uexact, 2) / норма (uexact, 2) = 3.1e-12

норма (U-uexact, инф) / норма (uexact, инф) = 2.7e-12

Казалось бы, решение подходит для 11-12 цифр, а номер условия порядка 1e8.

Однако ситуация с поэлементными ошибками более интересна.

макс (абс (U-uexact)) = 2.7e-12

Это все еще выглядит хорошо.

макс (абс ((и-uexact) ./ uexact) = 6.1e + 9

Ничего себе - есть очень большая относительная ошибка по крайней мере в одном компоненте решения.

Что произошло? Точное решение этого уравнения имеет крошечные компоненты (например, 1,9e-22), в то время как приблизительное решение достигает гораздо большего значения 9e-14. Это скрыто нормальным измерением относительной ошибки (будь то 2-норма или бесконечность-норма) и становится видимым только тогда, когда вы смотрите на элементарные относительные ошибки и берете максимум.

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


Как вы отметили в этом вопросе, число условий неособой матрицы дает оценку относительной ошибки в худшем случае для решения возмущенной системы уравнений. То есть, если мы точно решим A ( x + Δ x ) = b + Δ b и точно определим A x = b , тоκ(A)A(x+Δx)=б+ΔбAИксзнак равноб

Δxxκ(A)Δbb

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

Худшем случае ошибки случай имеет место , когда является левым особым вектором A , соответствующий наименьшей сингулярного значения A . Наилучший случай имеет место, когда ΔΔbAA является левым особым вектором A , соответствующий наибольшему значению особой A . Когда Δ b является случайным, то вы должны смотреть на проекции Δ b на все левые особые векторы A и соответствующие особые значения. В зависимости от спектра А , дела могут идти очень плохо или очень хорошо. ΔbAAΔbΔbAA

Рассмотрим две матрицы , обе с условным числом 2-нормы 1.0 × 10 10 . Первая матрица имеет сингулярные значения 1 , 1 × 10 - 10 , , 1 × 10 - 10A1.0×101011×10101×1010 . Вторая матрица имеет особые значения , 1 , , 1 , 1 × 10 - 10 . 1111×1010

В первом случае случайное возмущение вряд ли будет в направлении первого левого сингулярного вектора и, более вероятно, будет близко к одному из сингулярных векторов с сингулярным значением . Таким образом, относительное изменение в решении, вероятно, будет очень большим. Во втором случае практически любое возмущение будет близко по направлению к сингулярному вектору с сингулярным значением 1 , и относительное изменение в решении будет небольшим. 1×10101

PS (добавлено позже, когда я вернулся с занятий йогой ...)

Формула для решения имеет видAΔx=Δb

Δx=VΣ1UTΔb=i=1nUiTΔbσiVi

По теореме Пифагора,

Δx22=i=1n(UiTΔbσi)2

Если мы продолжим , то эта сумма максимальна , когда Δ б = U н и свести к минимуму , когда Δ б = U 1 .Δb2=1Δb=UnΔb=U1

В ситуации , рассматриваемой здесь, является результатом случайных ошибок округления, так что U T я Д Ь значений все должны быть примерно той же величина. Слагаемые с меньшими значениями σ я внесет большой вклад в ошибку, в то время как члены с большими значениями сг я не будет способствовать много. В зависимости от спектра, это может быть намного меньше, чем предел наихудшего случая. ΔbUiTΔbσiσi

Брайан Борхерс
источник
Не будет ли этот аргумент предполагает , что это возможно (даже если маловероятны) , чтобы достичь наихудшего случая , связанный из для матрицы в примере? AFAIU, основываясь на моем ответе и на основе документации, это не должно быть возможным. κ(A)?getrs
Кирилл
@BrianBorchers Не могли бы вы пояснить, почему «Ошибка наихудшего случая возникает, когда является левым единственным вектором A, соответствующим наименьшему сингулярному значению A. Наилучший случай возникает, когда Δ b является левым единственным вектором A, соответствующим наибольшему сингулярное значение A «. имеет место? Из приведенного ниже примера это логично, но мне понадобятся некоторые формулы. Пусть СВД из А будет = U Σ V T . В первом случае, = Δ б σ 1 v T 1 +ΔbAAΔbAAAA=UΣVTA=Δbσ1v1T+i=2NuiσiviT. How to proceed?
Zoltán Csáti
A matrix, but the general effect is similar- unless you get really unlucky in the round-off errors, you typically do somewhat better than the pessimistic worst-case bound.
Brian Borchers
(-1) The discussion of component-wise relative errors in the output is seriously misleading.
Kirill
1

Т.Л., д - р Они сообщили в число обусловленности, но не обязательно правильное число обусловленности для матрицы, потому что есть разница.

*getrs

xx0xcond(A,x)ucond(A)u.
Here cond(A,x) is not quite the usual condition number κ(A), but rather
cond(A,x)=|A1||A||x|x,cond(A)=|A1||A|.
(Here inside the norm these are component-wise absolute values.) See, for example, Iterative refinement for linear systems and LAPACK by Higham, or Higham's Accuracy and Stability of Numerical Algorithms (7.2).

For your example, I took a pseudospectral differential operator for a similar problem with n=128, and there is in fact a big difference between |A1||A| and κ(A), I computed 7×103 and 2.6×107, which is enough to explain the observation that this happens for all right hand sides, because the orders of magnitudes roughly match what is seen in Table 3.1 (3-4 orders better errors). This doesn't work when I try the same for just a random ill-conditioned matrix, so it has to be a property of A.

An explicit example for which the two condition numbers don't match, which I took from Higham (7.17, p.124), due to Kahan is

(2111ϵϵ1ϵϵ),(2+2ϵϵϵ).
Another example I found is just the plain Vandermonde matrix on [1:10] with random b. I went through MatrixDepot.jl and some other ill-conditioned matrices also produce this type of result, like triw and moler.

Essentially, what's going on is that when you analyze the stability of solving linear systems with respect to perturbations, you first have to specify which perturbations you are considering. When solving linear systems with LAPACK, this error bound considers component-wise perturbations in A, but no perturbation in b. So this is different from the usual κ(A)=A1A, which considers normwise perturbations in both A and b.

Consider (as a counterexample) also what would happen if you don't make the distinction. We know that using iterative refinement with double precision (see link above) we can get the best possible forward relative error of O(u) for those matrices with κ(A)1/u. So if we consider the idea that linear systems can't be solved to accuracy better than κ(A)u, how would refining solutions possibly work?

P.S. It matters that ?getrs says the computed solution is the true solution of (A + E)x = b with a perturbation E in A, but no perturbation in b. Things would be different if perturbations were allowed in b.

Edit To show this working more directly, in code, that this is not a fluke or a matter of luck, but rather the (unusual) consequence of two condition numbers being very different for some specific matrices, i.e.,

cond(A,x)cond(A)κ(A).
function main2(m=128)
    A = matrixdepot("chebspec", m)^2
    A[1,:] = A[end,:] = 0
    A[1,1] = A[end,end] = 1
    best, worst = Inf, -Inf
    for k=1:2^5
        b = randn(m)
        x = A \ b
        x_exact = Float64.(big.(A) \ big.(b))
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        best, worst = min(best, err), max(worst, err)
    end
    @printf "Best relative error:       %.3e\n" best
    @printf "Worst relative error:      %.3e\n" worst
    @printf "Predicted error κ(A)*ε:    %.3e\n" cond(A, Inf)*eps()
    @printf "Predicted error cond(A)*ε: %.3e\n" norm(abs.(inv(A))*abs.(A), Inf)*eps()
end

julia> main2()
Best relative error:       2.156e-14
Worst relative error:      2.414e-12
Predicted error κ(A)*ε:    8.780e-09
Predicted error cond(A)*ε: 2.482e-12

Edit 2 Here is another example of the same phenomenon where the different conditions numbers unexpectedly differ by a lot. This time,

cond(A,x)cond(A)κ(A).
Here A is the 10×10 Vandermonde matrix on 1:10, and when x is chosen randomly, cond(A,x) is noticably smaller than κ(A), and the worst case x is given by xi=ia for some a.
function main4(m=10)
    A = matrixdepot("vand", m)
    lu = lufact(A)
    lu_big = lufact(big.(A))
    AA = abs.(inv(A))*abs.(A)
    for k=1:12
        # b = randn(m) # good case
        b = (1:m).^(k-1) # worst case
        x, x_exact = lu \ b, lu_big \ big.(b)
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        predicted = norm(AA*abs.(x), Inf)/norm(x, Inf)*eps()
        @printf "relative error[%2d]    = %.3e (predicted cond(A,x)*ε = %.3e)\n" k err predicted
    end
    @printf "predicted κ(A)*ε      = %.3e\n" cond(A)*eps()
    @printf "predicted cond(A)*ε   = %.3e\n" norm(AA, Inf)*eps()
end

Average case (almost 9 orders of magnitude better error):

julia> T.main4()
relative error[1]     = 6.690e-11 (predicted cond(A,x)*ε = 2.213e-10)
relative error[2]     = 6.202e-11 (predicted cond(A,x)*ε = 2.081e-10)
relative error[3]     = 2.975e-11 (predicted cond(A,x)*ε = 1.113e-10)
relative error[4]     = 1.245e-11 (predicted cond(A,x)*ε = 6.126e-11)
relative error[5]     = 4.820e-12 (predicted cond(A,x)*ε = 3.489e-11)
relative error[6]     = 1.537e-12 (predicted cond(A,x)*ε = 1.729e-11)
relative error[7]     = 4.885e-13 (predicted cond(A,x)*ε = 8.696e-12)
relative error[8]     = 1.565e-13 (predicted cond(A,x)*ε = 4.446e-12)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

Worst case (a=1,,12):

julia> T.main4()
relative error[ 1]    = 0.000e+00 (predicted cond(A,x)*ε = 6.608e-13)
relative error[ 2]    = 1.265e-13 (predicted cond(A,x)*ε = 3.382e-12)
relative error[ 3]    = 5.647e-13 (predicted cond(A,x)*ε = 1.887e-11)
relative error[ 4]    = 8.895e-74 (predicted cond(A,x)*ε = 1.127e-10)
relative error[ 5]    = 4.199e-10 (predicted cond(A,x)*ε = 7.111e-10)
relative error[ 6]    = 7.815e-10 (predicted cond(A,x)*ε = 4.703e-09)
relative error[ 7]    = 8.358e-09 (predicted cond(A,x)*ε = 3.239e-08)
relative error[ 8]    = 1.174e-07 (predicted cond(A,x)*ε = 2.310e-07)
relative error[ 9]    = 3.083e-06 (predicted cond(A,x)*ε = 1.700e-06)
relative error[10]    = 1.287e-05 (predicted cond(A,x)*ε = 1.286e-05)
relative error[11]    = 3.760e-10 (predicted cond(A,x)*ε = 1.580e-09)
relative error[12]    = 3.903e-10 (predicted cond(A,x)*ε = 1.406e-09)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

Edit 3 Another example is the Forsythe matrix, which is a perturbed Jordan block of any size of the form

A=(010000100001ϵ000).
This has A=1, A1=ϵ1, so κ(A)=ϵ1, but |A1|=A1=|A|1, so cond(A)=1. And as can be verified by hand, solving systems of linear equations like Ax=b with pivoting is extremely accurate, despite the potentially unbounded κ(A). So this matrix too will yield unexpectedly precise solutions.

Edit 4 Kahan matrices are also like this, with cond(A)κ(A):

A = matrixdepot("kahan", 48)
κ, c = cond(A, Inf), norm(abs.(inv(A))*abs.(A), Inf)
@printf "κ=%.3e c=%.3e ratio=%g\n" κ c (c/κ)

κ=8.504e+08 c=4.099e+06 ratio=0.00482027
Kirill
источник
The condition numbers in the paper referred to by the OP are two-norm condition numbers. If you go back to reference [17] by ElBarbary you'll see that in the earlier paper these were two-norm condition numbers. Also, I setup the examples from this paper using DMsuite, and got nearly exactly the same 2-norm condition numbers as reported in the paper.
Brian Borchers
The infinity norm condition norm numbers for these examples that I got using dmsuite and Chebyshev interpolation were similar in magnitude to the two-norm condition numbers. I don't think that the difference between 2-norm in infinity-norm condition numbers is that important for this particular example.
Brian Borchers
I believe that the errors reported in the paper are absolute rather than relative errors (it doesn't make too much difference except for ϵ=0.01, where the solution dips down close to 0.
Brian Borchers
For ϵ=0.01 and N=1024, the relative errors for the parts of the solution that are near 0 are huge, but the absolute errors are small. I agree that the paper was very vague about what condition number was used and about what the "errors" were exactly (relative or absolute errors.)
Brian Borchers
@BrianBorchers I'm not sure what you mean: this isn't the difference between 2-norm and infty-norm condition numbers, but rather normwise- and component-wise condition numbers (component-wise relative perturbations in the input, not component-wise relative errors in the output as in your answer).
Kirill