【问题标题】:Algorithm for finding Eigenvectors given Eigenvalues of a 3x3 matrix in C#在C#中给定3x3矩阵的特征值查找特征向量的算法
【发布时间】:2018-11-25 17:22:22
【问题描述】:

我正在尝试使用 PCA 为我的网格找到最佳的 OOBB hitbox。为了做到这一点,我需要特征向量,但我有点迷失如何在不使用庞大库的情况下计算它们。

我实现了一种算法,该算法在给定 3x3 矩阵的情况下计算三个特征值。这个代码最初来自Wikipedia

private Vector3 CalculateEigenvalues(ref Matrix3 A)
{
    Vector3 val = new Vector3(0, 0, 0);

    float p1 = A.M12 * A.M12 + A.M13 * A.M13 + A.M23 * A.M23;
    if (p1 == 0)
    {
        val.X = A.M11;
        val.Y = A.M22;
        val.Z = A.M33;
    }
    else
    {
        float q = A.Trace / 3f;
        float p2 = (float)(Math.Pow(A.M11 - q, 2) + Math.Pow(A.M22 - q, 2) + Math.Pow(A.M33 - q, 2)) + 2 * p1;
        float p = (float)Math.Sqrt(p2 / 6);
        Matrix4 I = Matrix4.Identity;
        Matrix4.Mult(ref I, q, out Matrix4 tmp);
        Matrix4 tmp2 = Matrix4.Subtract(new Matrix4(A), tmp);
        Matrix4 B = Matrix4.Mult(tmp2, 1 / p);
        float r = new Matrix3(B).Determinant / 2;

        float phi = 0;
        if (r <= -1)
            phi = (float)Math.PI / 3;
        else if (r >= 1)
            phi = 0;
        else
            phi = (float)Math.Acos(r) / 3;
        val.X = q + 2 * p * (float)Math.Cos(phi);
        val.Z = q + 2 * p * (float)Math.Cos(phi + (2 * Math.PI / 3));
        val.Y = 3 * q - val.X - val.Z;
    }
    return val;
}

但是,Wikipedia 文章没有用于计算三个特征值的特征向量的代码。我试图理解这个话题,但我的数学技能非常有限。我必须在每个教程中搜索每个单词。

所以我的问题是:

如果我有 3x3 矩阵和三个特征值,是否有任何简单的方法可以在不使用外部库的情况下计算相应的特征向量?

【问题讨论】:

  • 您需要求解一个由三个线性方程组成的系统。这本身就是一项非常简单的任务,但需要一些编码。
  • 好吧,我想我可以想出一个算法来解决一个三线性方程组。但是我怎么知道要放入哪些方程呢?
  • 对于每个特征值 λ,对应的特征向量 vAv = λv 的解。这扩展为 3 个方程——向量的每个分量一个。

标签: c# algorithm eigenvalue eigenvector


【解决方案1】:

当所有特征值都具有相同的algebraic and geometric 多重性(旋转矩阵就是这种情况)时,绝对简单的实现是这样的:

// Observe that the function doesn't use rZ,
// it is expected that it will become zero vector in triangular form
static Vector3 EigenVector(Vector3 rX, Vector3 rY, Vector3 rZ, float lambda)
{
    // Move RHS to LHS
    rX.X -= lambda;
    rY.Y -= lambda;
    // Transform to upper triangle
    rY -= rX * (rY.X / rX.X);
    // Backsubstitute
    var res = new Vector3(1f);
    res.Y = -rY.Z / rY.Y;
    res.X = -(rX.Y * res.Y + rX.Z * res.Z) / rX.X;
    return res;
}

// Case of eigenvalue with algebraic multiplicity two
static (Vector3, Vector3) EigenVector2(Vector3 rX, Vector3 rY, Vector3 rZ, float lambda)
{
    // Move RHS to LHS
    rX.X -= lambda;
    float x2 = rX.Y / rX.X;
    float x3 = rX.Z / rX.X;
    return (new Vector3(x2, 1, 0), new Vector3(x3, 0, 1));
}

static void Main(string[] args)
{
    var rX = new Vector3(1, -3, 3);
    var rY = new Vector3(3, -5, 3);
    var rZ = new Vector3(6, -6, 4);
    var e = EigenVector(rX, rY, rZ, 4);
    var e2 = EigenVector2(rX, rY, rZ, 2);

    System.Diagnostics.Debug.WriteLine(e.ToString());
    System.Diagnostics.Debug.WriteLine(e2.Item1.ToString());
    System.Diagnostics.Debug.WriteLine(e2.Item2.ToString());
}



在现实生活中,需要进行大量的错误检查。输入数据取自this paper

【讨论】:

  • 那么 lamdba 是一个特征值,而 rX、rY 和 rZ 向量是矩阵的行?还是它的列?
  • @3DLearner 是的,lambda 是一个特征值和 rX, rY, rZ - 行,我在 C# 中没有找到 Matrix3 类。
  • 谢谢!我将对此进行测试。当我完成测试时,我会回复你。顺便说一句:Matrix3 类来自我在 C# 中用于 OpenGL 的 OpenTK 库。
  • 在您的实现中,您需要检查除以零。比如说,如果在第一个函数中 rX.X 非常接近于零,则需要交换行。像这样的算法很容易在数值上变得不稳定。如果您需要做更多的数学运算,最好使用经过良好测试的外部库。它将帮助您更快地构建应用程序,并且不会被血淋淋的细节所淹没。
  • 因为我要求最简单的方法来做到这一点,我很高兴将此问题标记为已解决。非常感谢您的帮助!我现在终于有了更好的 Hitboxes ;-)
猜你喜欢
  • 2022-11-24
  • 2018-10-31
  • 2017-10-21
  • 1970-01-01
  • 2012-08-08
  • 2019-04-27
  • 2013-05-15
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多