无法使用 OjAlgo 执行正定对称矩阵的广义特征值问题 - 我在做什么错?

Cannot perform generalized eigenvalue problem with positive definitive symmetrical matricies with OjAlgo - What wrong am I doing?

正在尝试求解以下形式的广义特征值:

A*V = B*V*D

通过使用 OjAlgo。根据文档 here AB 半身像是实对称或复厄密矩阵,B 是正定的。在这种情况下,AB 都是对称的和肯定的。

OjAlgo 是唯一可以解决广义特征值问题的 Java 数学库。所以这必须工作。但为什么我的输出说我无法解决它?

public class Eig {

    static Logger logger = LoggerFactory.getLogger(Eig.class);

    // A*V = B*D*V - Find D and V - Will not work for current OjAlgo version
    static public void eig(MatrixStore<Double> Sb, MatrixStore<Double> Sw, Primitive64Store D, Primitive64Store V, long dim) {

        // Create eigA and eigB from symmetrical positive definitive A and B
        Primitive64Matrix eigA = Primitive64Matrix.FACTORY.rows(Sb.toRawCopy2D());
        Primitive64Matrix eigB = Primitive64Matrix.FACTORY.rows(Sw.toRawCopy2D());

        System.out.println("Check if eigA and eigB are symmetrical:");
        System.out.println(eigA.isSymmetric());
        System.out.println(eigB.isSymmetric());

        System.out.println("Check if eigA and eigB are positive definitive:");
        Primitive64Matrix z = Primitive64Matrix.FACTORY.makeFilled(dim, 1, new Weibull(5, 2));
        System.out.println("Positive definitive:");
        System.out.println(z.transpose().multiply(eigA).multiply(z).get(0, 0)); // z^T*eigA*z
        System.out.println(z.transpose().multiply(eigB).multiply(z).get(0, 0)); // z^T*eigB*z

        // Perform [A][V] = [B][V][D]
        Eigenvalue.Generalised<Double> eig = Eigenvalue.PRIMITIVE.makeGeneralised(eigA, Generalisation.A_B);
        boolean success = eig.computeValuesOnly(eigA, eigB);
        if (success == false)
            logger.error("Could not perform eigenvalue decomposition!");


        System.out.println("Check if D and V are null");
        System.out.println(eig.getD() == null);
        System.out.println(eig.getV() == null);

        // Copy over to D, V
        D.fillColumn(0, eig.getD().sliceDiagonal());
        double[][] eigV = eig.getV().toRawCopy2D();
        for (int i = 0; i < dim; i++) {
            V.fillRow(i, Access1D.wrap(eigV[i]));
        }

        // Sort eigenvalues and eigenvectors descending by eigenvalue
        if(eig.isOrdered() == false)
            Sort.sortdescended(V, D, dim); 

    }
}

我错过了什么?

Check if eigA and eigB are symmetrical:
true
true
Check if eigA and eigB are positive definitive:
Positive definitive:
1.0766814686417156E10
1.1248634208301022E9
Check if D and V are null:
true
true

更新 1:

如果 AB 都只有正实数值,该过程将起作用。

        Primitive64Store mtrxA = Primitive64Store.FACTORY.makeSPD((int) dim);
        Primitive64Matrix eigA = Primitive64Matrix.FACTORY.rows(mtrxA.toRawCopy2D());
        Primitive64Store mtrxB = Primitive64Store.FACTORY.makeSPD((int) dim);
        Primitive64Matrix eigB = Primitive64Matrix.FACTORY.rows(mtrxB.toRawCopy2D());

        PrintMatrix.printMatrix(eigB);

        /*
         * There are several generalisations. 3 are supported by ojAlgo, specified by the enum:
         * Eigenvalue.Generalisation This factory method returns the most common alternative.
         */
        Eigenvalue.Generalised<Double> generalisedEvD = Eigenvalue.PRIMITIVE.makeGeneralised(eigA);
        // Generalisation: [A][V] = [B][V][D]

        // Use 2-args alternative

        generalisedEvD.decompose(eigA, eigB);

        System.out.println("Check if D and V are null");
        System.out.println(generalisedEvD.getD() == null); // false
        System.out.println(generalisedEvD.getV() == null); // false

更新二:

我用 GNU Octave 做了 运行 测试,似乎所有的特征值都是正的,其余的都是负的,但非常接近于零。

这是一个输出。它与我在 GNU Octave 和 OjAlgo 中使用的数据相同。我认为 e-18 可以算作零。

我构建了我的 AB,因为它们应该是对称的和确定的。这是由浮动值引起的吗?

   2.7414e+04
   9.4155e+03
   4.1295e+03
   3.1429e+03
  -8.4338e-16
  -1.6409e-15
          Inf
          Inf
          Inf
          Inf
          Inf
          Inf
          Inf
          Inf
          Inf
          Inf
   3.4910e-15
  -8.7739e-16
  -3.1775e-15
  -2.8213e-18
  -5.0274e-16
   1.7329e-18
  -1.1330e-15
   3.1024e-18
   2.3226e-15
  -1.6151e-16
  -6.8453e-16
   1.6111e-17
  -1.7850e-18
  -1.3411e-18
  -2.3916e-18

在您调用的第一个代码示例中:

eig.computeValuesOnly(eigA, eigB);

这只会给你特征值(没有向量或矩阵)。在第二个示例中,您改为调用通常的:

generalisedEvD.decompose(eigA, eigB);