Поиск корней с помощью сопутствующей матрицы

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

Я знаю, как превратить полином в сопутствующую матрицу, и я нашел скрипт, который выполняет QR-разложение: http://quantstart.com/articles/QR-Decomposition-with-Python-and-NumPy

И вот тут я потерялся: что делать дальше? Я думаю, что я должен рассчитать несколько разложений, но когда я делаю, я всегда получаю один и тот же результат (очевидно). Я также читал, что было бы полезно сначала преобразовать сопутствующую матрицу в форму Гессенберга - но как? Тогда есть "сдвиги" - что это?

Я также нашел http://www.nr.com/webnotes/nr3web17.pdf но так как я не понимаю ничего из этого, я хотел бы знать, есть ли какой-нибудь более простой метод (даже если он медленнее или менее стабилен).

Другими словами: чтение http://en.wikipedia.org/wiki/QR_algorithm

  • "пусть A будет реальной матрицей, из которой мы хотим вычислить собственные значения"
    хорошо, это моя сопутствующая матрица, верно?

  • "вычисляем QR-разложение Ak=QkRk"
    это было бы Q, R = householder(A) с самой первой ссылки, верно?

  • "Затем мы формируем Ak+1 = RkQk"
    легко, просто умножьте R и Q

  • "При определенных условиях [2] матрицы Ak сходятся к треугольной матрице, форма Шура для A. Собственные значения треугольной матрицы перечислены на диагонали, и задача на собственные значения решена".
    ...чего ждать? Я старался:

    for i in range(100):
    Q, R = householder(A)
    A = mult_matrix(R, Q)

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

Пожалуйста, кто-нибудь может мне это объяснить?

PS: я не хочу слепо использовать LAPACK или подобное, так как хочу понять, как это работает, по крайней мере, в очень упрощенном виде.

PPS: есть также http://adorio-research.org/wordpress/?p=184 (хотя не уверен, чем он отличается от первого метода...)

Источник

1 ответ

Решение

Если ваша сопутствующая матрица является матрицей преобразования коэффициента линейной полиномиальной операции, которая отображает q(x) в x*q(x) mod p(x)

где p(x)=x^(n+1)+p_n*x^n+...+p_1*x+p_0,

Ясно, что А имеет форму

0 0 0 ... 0 -p_0
1 0 0 ... 0 -p_1
0 1 0 ... 0 -p_2
. . . ... . . . 
0 0 0 ... 1 -p_n

который уже находится в форме Гессенберга. Поскольку эта форма сохраняется во время алгоритма QR, вы можете использовать декомпозицию QR с вращениями Гивенса, когда происходят только вращения, близкие к диагонали.

В несмещенном алгоритме QR вы должны наблюдать заметное развитие, по крайней мере, в правом нижнем блоке 3x3.

Если вы берете одно из собственных значений нижнего правого блока 2x2 и вычитаете его из каждого диагонального элемента, то вы получаете смещенные QR-алгоритмы в соответствии с Фрэнсисом. Сдвиг s, вычитаемое число, является текущей наилучшей оценкой наименьшего собственного значения. Обратите внимание, что вполне вероятно, что теперь вы покинули реальный домен и должны выполнять вычисления со сложными матричными записями. Вы должны сохранить сдвиг в памяти и добавить к нему любой новый сдвиг в последующих шагах и добавить комбинированный сдвиг обратно к любому найденному собственному значению.

Разбиение матрицы происходит всякий раз, когда любой из поддиагональных элементов практически равен нулю. Если разделение происходит в последнем ряду, то последний диагональный элемент является собственным значением сдвинутой матрицы (As*I). Если разделение разделяет последний блок 2x2, то можно легко определить его собственные значения, которые снова являются собственными значениями сдвинутой матрицы.

Если разделение происходит где-либо еще, алгоритм QR рекурсивно применяется к диагональным блокам отдельно.

  • Приложение I

Сходимость всех вариантов алгоритма измеряется записями под диагональю, сходящимися к нулю.

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

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

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

Приложение II

Это линия

   alpha = -cmp(x[0],0) * norm(x)

в процедуре домохозяина. В большинстве ситуаций x[0] не будет точно 0, так что это производит знак. Однако в случае сопутствующей матрицы x[0]==0 по построению, так что никаких признаков не производится, альфа получает неправильное значение 0. Измените его на более пешеходный

    alpha = norm(x)
    if x[0] < 0: alpha = -alpha

и это работает хорошо.

def companion_matrix(p):
    n=len(p)-1
    C=[[float(i+1 == j) for i in xrange(n-1)] for j in xrange(n)]
    for j in range(n): C[j].append(-p[n-j]/p[0])
    return C


def QR_roots(p):
    n=len(p)-1
    A=companion_matrix(p)
    for k in range(10+n):
        print "step: ",k+1," after ",(k+1)*(5+n), "iterations"
        for j in range(5+n):
            Q,R=householder(A)
            A=mult_matrix(R,Q)
        print("below diagonal")
        pprint([ A[i+1][i] for i in range(n-1) ])
        print("diagonal")
        pprint([ A[i][i] for i in range(n) ])
        print("above diagonal")
        pprint([ A[i][i+1] for i in range(n-1) ])


p=[ 1, 2, 5, 3, 6, 8, 6, 4, 3, 2, 7]
QR_roots(p)

#for a case with multiple roots at 1 and 4
#p= [1,-11,43,-73,56,-16]
#QR_roots(p)
Источник
Другие вопросы по тегам