21

Eigen是否具有存储密集、固定大小、对称矩阵的有效类型?(嘿,它们无处不在!)

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

如果没有这样的事情,我应该采取哪些行动(看起来像这样)将这种类型引入 Eigen?它有“视图”的概念吗?我可以为我自己的类型写一些像“矩阵视图”这样的东西,这会使它成为 Eigen-friednly 吗?

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

4

3 回答 3

12

对称矩阵的打包存储是矢量化代码的一敌人,即速度。标准做法是将相关的 N*(N+1)/2 系数存储在全密集 NxN 矩阵的上三角或下三角部分中,而剩余的 (N-1)*N/2 不被引用。然后通过考虑这种特殊的存储来定义对称矩阵上的所有操作。在 eigen 中,您有三角和自伴视图的概念来获得这一点。

特征参考:(对于实矩阵selfadjoint ==对称)。

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

除非内存是一个大问题,否则我建议将矩阵的未引用部分留空。(更易读的代码,没有性能问题。)

于 2012-11-15T21:39:47.523 回答
11

是的,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;

不过它很麻烦,如果你的记忆真的很珍贵,我只会使用它。

于 2012-11-15T20:48:09.287 回答
1

对称矩阵的高效类型

您只需将值分配给矩阵的下/上三角部分并使用本征三角和自伴随视图。但是,我已经在小型固定大小的矩阵上进行了测试。我注意到在性能方面,使用视图并不总是最好的选择。考虑以下代码:

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矩阵来说甚至更快)。

一些替代品。

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)只存储三角形部分并实现自己的操作

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

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产品见下例(性能与方案一中的产品相同)。该算法只访问矩阵的下三角部分。

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;
}

Matrix2d*Vector2d在我的例子中,与使用完整矩阵(2x2 = 4 个元素)的产品相比,产品的性能增益为10%。

于 2019-02-23T09:45:06.110 回答