1

一変数多項式の実根をすべて見つけたいと思います。たとえば、Jenkins-Traub アルゴリズムを使用できますが、コンパニオン マトリックスを使用してそれを解決する方法を学びたいと考えています。

多項式をコンパニオン行列に変換する方法を知っており、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 の Schur 形式である三角行列に収束します。三角行列の固有値は対角線上にリストされ、固有値の問題が解決されます。」
    ...待って、何？私は試した：

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

しかし、何の進歩も見られず、正しい根に近い数字さえ見当たりません。

お願いします、誰か私にこれを説明してもらえますか?

PS: 少なくとも非常に単純化された用語で、LAPACK の仕組みを理解したいので、やみくもに LAPACK などを使用したくありません。

PPS: http://adorio-research.org/wordpress/?p=184もあります(最初の方法とどう違うのかわかりませんが...)

4

1 に答える 1

1

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 ブロックの固有値の 1 つを取り、それをすべての対角要素から差し引くと、Francis によるシフトされた QR アルゴリズムが得られます。減算された数であるシフト s は、最小の固有値の現在の最良の推定値です。おそらく、実際のドメインを離れて、複雑な行列エントリを計算する必要があることに注意してください。シフト s をメモリに保持し、後続のステップで新しいシフトを追加し、組み合わせたシフトを見つかった固有値に追加する必要があります。

サブ対角エントリのいずれかが実質的にゼロになると、行列の分割が発生します。分割が最後の行で発生した場合、最後の対角要素はシフトされた行列の固有値 (As*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構造上、符号が生成されないため、alpha は間違った値の 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)
于 2013-12-11T09:22:28.820 に答える