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

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

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

В явном виде A имеет вид

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, вычитаемое число, является текущей наилучшей оценкой наименьшего собственного значения. Обратите внимание, что, вполне вероятно, вы покинули реальную область и должны выполнять вычисления со сложными матричными элементами. Вы должны сохранить сдвиг в памяти и добавить к нему любой новый сдвиг на последующих шагах, а также добавить комбинированный сдвиг обратно к любому найденному собственному значению.

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

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

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

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

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

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

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

Приложение 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)