将shift和Hessenberg实现为一个已经运行(缓慢)的QR算法

Question

TLDR：帮助我让底部的函数做我认为它应该做的事情。

因此，在一次面试不顺利后，我决定写自己的库，这样，如果我再次被问及矩阵乘法，我会做好更好的准备。

我知道用Java做这件事会牺牲一些性能（至少我认为我是这样做的）。此外，我知道已经有一些图书馆是为这个目的而建的。关键是要学会足够的知识，能够即兴处理我可能需要的任何矩阵算法，并学会在哪里找到我不需要的答案。我用Java做这件事是为了充实我对这门语言的理解。

除了最后一部分：特征值之外，我已经完成了在库中要做的大部分工作。

我已经实现了一个相当基本的QR分解方法，它似乎很有效（阅读：将我库中的一堆随机矩阵与可信计算器的输出进行比较）。

问题是它的数量级太慢了。对于一个50x50的矩阵，几乎需要一分钟才能得到特征值。

现在，这是因为我已经将迭代次数设置为500次，以处理我在测试时遇到的一个病理病例。我想将迭代次数减少到10次左右，这几乎足以使每个矩阵收敛。

我浏览了一下网页，发现了一本关于高级二维码方法的书。基本上，根据我对这篇论文的理解，如果你首先将矩阵转换为上海森堡形式，它应该会更快地收敛几个数量级。此外，如果你实现了移位，它应该是二次收敛的。

我阅读并实现了这两种算法。它们确实使收敛更快；问题是，现在它没有通过上一个问题测试用例，默默地吐出了错误的特征值。

编辑：我测试了Hessenberg函数，它可以处理任意矩阵；在中，它输出的数字与MATLAB和wolfram相同。然而，当我将其作为第一步添加时，迭代次数会增加。我现在正在数学堆栈交换上询问底层算法，但转换确实是让我头疼的部分

我不明白为什么。通过我有限的数学理解，这本书说，如果矩阵a经过Hessenberg变换H（a），它应该仍然具有相同的特征值。换挡也是如此。但当我实现其中一个或两个算法时，特征值会发生变化。

我的问题是要么我错误地实现了算法，要么我误解了它背后的数学

我所说的论文供参考：http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter3.pdf

编辑：存储库链接（以及此代码所依赖的其余代码）：https://github.com/rwthompsonii/matrix-java

编辑：按照规则的相关功能：

public static Complex[] eigenvalues(SquareMatrix A) {
    Complex[] e = new Complex[A.getRows()];
    QRDecomposition qr = new QRDecomposition();
    qr.iterations = 0;
    int total_iter = 0;
    int num_eigen_found = 0;
    SquareMatrix QRIterator = new SquareMatrix(A);
    //in general, QR decomposition will converge faster from an upper
    //Hessenberg matrix.  so, first things first, we bring QRIterator to that form
    //QRIterator = new SquareMatrix(qr.hessenberg(QRIterator));
    int max = MAX_ITERATIONS;
    //double lastElement;
    //SquareMatrix ScaledIdentity;
    do {
        System.out.println("Pre-decompose: QRIterator (Iteration#" + (qr.iterations + 1) + "):n" + QRIterator);
        if (QRIterator.getRows() == 1) {
            //very last 1x1 element in matrix
            e[num_eigen_found++] = new Complex(
                    QRIterator.getMatrix()[0][0]
            );
            break;
        } else {
            /*lastElement = QRIterator.getMatrix()[QRIterator.getRows() - 1][QRIterator.getColumns() - 1];
            ScaledIdentity = new SquareMatrix(Matrix.IdentityMatrix(QRIterator.getRows()).scale(lastElement));
            try {
                QRIterator = new SquareMatrix(QRIterator.subtract(ScaledIdentity));
            } catch (DimensionMismatchException ex) {
                System.out.println("Unexpected execption during QRIterator -= I*alpha, bailing.");
                System.exit(-1);
            }*/
            qr.decompose(QRIterator);
        }
        try {
            QRIterator = new SquareMatrix(qr.R.mult(qr.Q)/*.add(ScaledIdentity)*/);
        } catch (DimensionMismatchException ex) {
            System.out.println("An unexpected exception occurred during QRIterator = R*Q, bailing.");
            System.exit(-1);
        }
        qr.iterations++;
        //testing indicates that MAX_ITERATIONS iterations should be more than sufficient to converge, if its going to at all
        if (qr.iterations == max || Math.abs(QRIterator.getMatrix()[QRIterator.getRows() - 1][QRIterator.getColumns() - 2]) < CONVERGENCE_CHECK) {
            System.out.println("QRIterator (at max iteration or converged) (Iteration#" + (qr.iterations + 1) + "):n" + QRIterator);
            if (Math.abs(QRIterator.getMatrix()[QRIterator.getRows() - 1][QRIterator.getColumns() - 2]) < CONVERGENCE_CHECK) {
                //then the value at M[n][n] is an eigenvalue and it is real
                e[num_eigen_found++] = new Complex(
                        QRIterator.getMatrix()[QRIterator.getRows() - 1][QRIterator.getColumns() - 1]
                );
                //System.out.println("e[" + (num_eigen_found - 1) + "]:t" + e[num_eigen_found - 1] + "nQRIterator before deflation:n" + QRIterator);
                double[][] deflatedMatrix = deflate(QRIterator.getMatrix(), 1);
                QRIterator = new SquareMatrix(deflatedMatrix);
                //System.out.println("nQRIterator after deflation:n" + QRIterator);
                total_iter += qr.iterations;
                qr.iterations = 0;  //reset the iterations counter to find the next eigenvalue
            } else {
                //this is a 2x2 matrix with either real or complex roots.  need to find them.
                //characteristic equation of 2x2 array => E^2 - (w + z)E + (wz - xy) = 0 where E = eigenvalue (possibly pair, possibly singular, possibly real, possibly complex)
                // and the matrix {{w, x}, {y, z}} is the input array, the task is to calculate the root(s) of that equation
                //that is a quadratic equation => (root = (-b +- sqrt(b^2  - 4ac))/2a)
                //determinant b^2 - 4ac will determine behavior of roots => positive means 2 real roots, 0 means 1 repeated real root, negative means conjugate pair of imaginary roots
                //first, get the wxyz from the (possibly bigger) matrix
                int n = QRIterator.getRows();
                double w = QRIterator.getMatrix()[n - 2][n - 2];
                double x = QRIterator.getMatrix()[n - 2][n - 1];
                double y = QRIterator.getMatrix()[n - 1][n - 2];
                double z = QRIterator.getMatrix()[n - 1][n - 1];
                //a not used since it's = 1
                double b = -(w + z);
                double c = (w * z - x * y);
                //calculate determinant of quadratic equation
                double determ = b * b - 4 * c;
                if (determ >= 0) {
                    //one or two real roots 
                    double sqrt_determ_real = Math.sqrt(determ);
                    e[num_eigen_found++] = new Complex((-b + sqrt_determ_real) / 2.0);
                    e[num_eigen_found++] = new Complex((-b - sqrt_determ_real) / 2.0);
                    //in the zero determinant case that's simply going to add the same eigenvalue to the list twice.  I'm ok with that for now.
                } else if (determ < 0) {
                    //conjugate pair of complex roots
                    double sqrt_determ_imag = Math.sqrt(-determ);
                    e[num_eigen_found++] = new Complex(-b / 2.0, sqrt_determ_imag / 2.0);
                    e[num_eigen_found++] = new Complex(-b / 2.0, -sqrt_determ_imag / 2.0);
                }
                if (QRIterator.getRows() > 2) {
                    total_iter += qr.iterations;
                    qr.iterations = 0;  //reset the iterations counter to find the next eigenvalue
                    double[][] deflatedMatrix = deflate(QRIterator.getMatrix(), 2);
                    QRIterator = new SquareMatrix(deflatedMatrix);
                } 
            }
        }
        //QRIterator = new SquareMatrix(qr.hessenberg(QRIterator));
    } while (qr.iterations < max);
    //used for debugging here
    /*System.out.println("Finished iterating.  Iterations:t" + total_iter
     + "nFinal value of qr.Q:n" + qr.Q + "nFinal value of qr.R:n" + qr.R
     + "nFinal value of QRIterator:n" + QRIterator
     + "nOriginal SquareMatrix A:n" + A);
     */
    return e;
}

编辑：在讽刺之前，我需要清理很多垃圾，比如我用来调试的一堆打印语句，主要是因为我不想通过500次迭代来查看值。我想让它工作起来，然后把它清理干净，以达到我自己相当不错的可读性标准。我知道需要进行一些重构，因为函数太长了。但首先，它需要工作。帮一个人？

Answer 1

我认为一个好的第一步是废除瑞利商移位法、威尔金森移位法和双移位法，直到你可以直接将其简化为Hessenburg形式，即你的第一个参考文献的算法3.3，第61页。

在这种情况下，迭代次数是固定的，一开始就有点简单。根据第59页的定义，计算Householder反射器P_k，其中k从1到n-2，包括向量范数和rho。请注意，在实际情况下，通常设置rho = -sign(x_1)，请参阅文本。

在纸上写一个简单的2x2示例（如果你感兴趣的话，可以写3x3），然后一步一步地print你的Java计算，以确保它们与你的手写作品一致。此外，请检查您是否遇到浮点问题等。

最后，如果您的代码没有轻而易举地找到任意50x50矩阵的特征值，请不要感到沮丧。这些缩减的大多数实现都经过了大量优化，很少会像文献中那样逐字翻译伪代码。如果你喜欢Python，可以看看其他库，比如Sage或NumPy；或者可能是C++的Boost。对于这种问题，使用这些可能会更有趣。