密集对称矩阵的特征有效型

Question

Eigen是否具有存储密集，固定大小，对称矩阵的有效类型？ （嘿，它们无处不在！）

即对于N = 9，它应该只存储（1 + 9）* 9/2 == 45个元素并且它具有适当的操作。例如，应该有效地添加两个对称矩阵，其返回相似的对称矩阵。

如果没有这样的东西，我应该做哪些动作（看起来像this）给Eigen？它有“观点”的概念吗？我可以为自己的类型写一些类似“矩阵视图”的东西，这会让它成为特征性的吗？

附：可能我可以使用map将普通数组视为1xN矩阵，并对其进行操作。但这不是最干净的解决方案。

Answer 1

是的，eigen3有views的概念。但它对存储没有任何作用。尽管如此，您可能能够为相同类型的两个对称矩阵共享更大的块：

Matrix<float,4,4> A1, A2; // assume A1 and A2 to be symmetric
Matrix<float,5,4> A;
A.topRightCorner<4,4>().triangularView<Upper>() = A1;
A.bottomLeftCorner<4,4>().triangularView<Lower>() = A2;

虽然它非常麻烦，但如果你的记忆非常珍贵，我只会使用它。

Answer 2

对称矩阵的Packed storage是矢量化代码的大敌，即速度。标准做法是将相关的N *（N + 1）/ 2系数存储在全密集NxN矩阵的上部或下部三角形部分中，并使剩余的（N-1）* N / 2不参考。然后通过考虑这种特殊的存储来定义对称矩阵上的所有操作。在本征中你有triangular and self-adjoint views的概念来获得这个。

来自eigen参考:(对于真实矩阵selfadjoint == symmetric）。

就像三角矩阵一样，您可以参考方阵的任何三角形部分，将其视为自相关矩阵，并执行特殊和优化的操作。同样，相反的三角形部分从未被引用，并且可用于存储其他信息。

除非记忆是一个大问题，否则我建议将矩阵的未引用部分留空。 （更易读的代码，没有性能问题。）

Answer 3

对称矩阵的有效类型

您只需将值分配给矩阵的下/上三角形部分，并使用特征三角和自相连视图。但是，我已经在小型固定尺寸矩阵上进行了测试。我注意到性能方​​面，使用视图并不总是最好的选择。请考虑以下代码：

Eigen::Matrix2d m;
m(0,0) = 2.0;
m(1,0) = 1.0;
// m(0,1) = 1.0;
m(1,1) = 2.0;
Eigen::Vector2d v;
v << 1.0,1.0;
auto result = m.selfadjointView<Eigen::Lower>()*v;

与下面提供的替代解决方案相比，最后一行中的产品相当慢（在我的情况下double 2x2矩阵慢了约20％）。 （使用完整矩阵的产品，通过取消注释m(0,1) = 1.0;，并使用auto result = m*v，对于double 2x2矩阵来说甚至更快）。

编辑：我忘记了别名。 auto result.noalise()加快了速度（下面的替代方案仍然更快）。

一些替代品。

1）在向量中存储对称矩阵

您可以将矩阵存储在大小为45的向量中。以向量格式对2个矩阵求和是很简单的（只需对向量求和）。但是你必须为产品编写自己的实现。

这是一个matrix * vector产品（密集，固定大小）的实现，其中矩阵的下半部分按列存储在向量中：

template <typename T, size_t S>
Eigen::Matrix<T,S,1> matrixVectorTimesVector(const Eigen::Matrix<T,S*(S+1)/2,1>& m, const Eigen::Matrix<T,S,1>& v)
{
    Eigen::Matrix<T,S,1> ret(Eigen::Matrix<T,S,1>::Zero());
    int counter(0);
    for (int i=0; i<S; ++i)
    {
        ret[i] += m(counter++)*v(i);
        for (int j=i+1; j<S; ++j)
        {
            ret[i] += m(counter)*v(j);
            ret[j] += m(counter++)*v(i);
        }
    }
    return ret;
}

2）仅存储三角形零件并执行您自己的操作

您当然也可以实现自己的产品matrix * vector，其中矩阵仅存储45个元素（假设我们存储下三角形部分）。这可能是最优雅的解决方案，因为它保持矩阵的格式（而不是使用表示矩阵的向量）。然后，您还可以使用以下示例中的Eigen函数：

template <typename T, size_t S>
Eigen::Matrix<T,S,S> symmMatrixPlusSymmMatrix( Eigen::Matrix<T,S,S>& m1, const Eigen::Matrix<T,S,S>& m2)
{
    Eigen::Matrix<T,S,S> ret;
    ret.template triangularView<Eigen::Lower>() = m1 + m2; // no performance gap here!
    return ret;
}

在上述函数中（2个对称矩阵的和），仅访问m1和m2的下三角形部分。请注意，triangularView在这种情况下不会出现性能差距（我根据我的基准确认这一点）。

至于matrix * vector产品，请参见下面的示例（与替代品1中的产品相同的性能））。该算法仅访问矩阵的下三角部分。

template <typename T, size_t S>
Eigen::Matrix<T,S,1> symmMatrixTimesVector(const Eigen::Matrix<T,S,S>& m, const Eigen::Matrix<T,S,1>& v)
{
    Eigen::Matrix<T,S,1> ret(Eigen::Matrix<T,S,1>::Zero());
    int counter(0);

    for (int c=0; c<S; ++c)
    {
        ret(c) += m(c,c)*v(c);
        for (int r=c+1; r<S; ++r)
        {
            ret(c) += m(r,c)*v(r);
            ret(r) += m(r,c)*v(c);
        }
    }
    return ret;
}

与使用全矩阵（2x2 = 4个元素）的产品相比，产品Matrix2d*Vector2d的性能提升在我的情况下为10％。