如何用Eigen将矩阵上三角视图映射为特征向量？

我正尝试根据一篇论文实现一款新型水库神经网络（类似循环神经网络，但仅训练输出层权重）。流程中需要从线性特征向量$O_{lin}$构造非线性特征向量$O_{nonlin}$，步骤如下：

• 二次多项式的所有单项式可通过外积$O_{lin} \otimes O_{lin}$获取
• 非线性特征向量由外积张量的上三角元素中的唯一单项式组成，论文用符号$\lceil\otimes\rceil$表示该操作。

举个例子：
线性特征向量：

1
2
3
4

对应的外积矩阵：

1  2  3  4
2  4  6  8
3  6  9 12
4  8 12 16

提取上三角元素后得到的非线性特征向量：

1
2
3
4
4
6
8
9
12
16

我尝试用Eigen实现该操作，编写的代码如下：

// 需要Eigen库
#include "Eigen/Dense"

#include <iostream>
#include <format>

int main() {
    Eigen::Vector<double, 4> o_lin;
    o_lin << 1, 2, 3, 4;

    const auto outer_prod{ o_lin * o_lin.transpose() };
    const auto upper_tri{ outer_prod.template triangularView<Eigen::Upper>() };
    Eigen::Vector<double, 16> upper_tri_vec;

    std::cout << outer_prod << "\n\n";
    std::cout << std::format( "upper_tri ({}, {}), size: {}\n",
                              upper_tri.rows(), upper_tri.cols(),
                              upper_tri.size() );
    std::cout << std::format( "upper_tri_vec ({}, {}), size: {}\n",
                              upper_tri_vec.rows(), upper_tri_vec.cols(),
                              upper_tri_vec.size() );

    // 运行时失败
    upper_tri.evalTo( upper_tri_vec );
    std::cout << upper_tri_vec << std::endl;

    // 无法编译，triangularView没有reshaped()方法
    // Eigen::Vector<double, 10> o_nonlin{ upper_tri.reshaped() };
}

运行时报错：

1  2  3  4
2  4  6  8
3  6  9 12
4  8 12 16

upper_tri (4, 4), size: 16
upper_tri_vec (16, 1), size: 16
Assertion failed: ((!(RowsAtCompileTime!=Dynamic) || (rows==RowsAtCompileTime)) && (!(ColsAtCompileTime!=Dynamic) || (cols==ColsAtCompileTime)) && (!(Rows
AtCompileTime==Dynamic && MaxRowsAtCompileTime!=Dynamic) || (rows<=MaxRowsAtCompileTime)) && (!(ColsAtCompileTime==Dynamic && MaxColsAtCompileTime!=Dynami
c) || (cols<=MaxColsAtCompileTime)) && rows>=0 && cols>=0 && "Invalid sizes when resizing a matrix or array."), function resize, file PlainObjectBase.h, l
ine 273.
zsh: abort      ./main

之后我自行实现映射操作，编写了模板函数：

// 简化类型名，std::ptrdiff_t == Eigen::Index
using Index = std::ptrdiff_t;

template <Weight T, Index R, Index C, Index N>
constexpr inline Eigen::Vector<T, N>
get_u_tri( const Eigen::Ref<Eigen::Matrix<T, R, C>> m ) {
    Vec<T, N> result;
    Index     pos{ 0 };
    for ( Index i{ 0 }; i < R; ++i ) {
        Index m_block_width{ std::max( C - i, Index{ 1 } ) };
        Index m_block_pos{ std::min( i, C - 1 ) };

        result.template segment( pos, m_block_width ) =
            m.template block( i, m_block_pos, 1, m_block_width );
        // result.template block( pos, 0, C - i, 1 ) =
        //     m.template block( i, m_block_pos, 1, m_block_width );

        pos += std::max( C - i, Index{ 1 } );
    }
    return result;
}

但仍出现运行时错误：

Assertion failed: (rows == this->rows() && cols == this->cols() && "DenseBase::resize() does not actually allow to resize."), function resize, file DenseB
ase.h, line 261.

请问如何解决这些问题？

错误原因分析 #

1. triangularView.evalTo()失败：triangularView是4x4的矩阵视图，而目标是16x1的向量，维度不匹配，Eigen的evalTo要求目标与视图维度完全一致，无法直接将二维视图转为一维向量。
2. 自行实现函数的问题：模板参数错误（Weight未定义）、Vec类型未定义、块宽度计算逻辑冗余、起始列判断多余，导致维度匹配失败触发断言。

正确实现方式 #

方式1：利用Eigen内置的对称矩阵打包视图（简洁高效）

外积矩阵是对称矩阵，可通过selfadjointView的packedView直接提取上三角元素的连续存储，再用Map转为向量：

#include "Eigen/Dense"
#include <iostream>

int main() {
    Eigen::Vector<double, 4> o_lin;
    o_lin << 1, 2, 3, 4;

    auto outer_prod = o_lin * o_lin.transpose();
    // 提取上三角打包视图，按列优先存储上三角元素
    auto packed_upper = outer_prod.selfadjointView<Eigen::Upper>().packedView();
    // 将打包后的内存映射为10维向量（4*(4+1)/2=10）
    Eigen::Vector<double, 10> o_nonlin = Eigen::Map<Eigen::Vector<double, 10>>(packed_upper.data());

    std::cout << o_nonlin << std::endl;
}

输出结果完全符合需求：

1
2
3
4
4
6
8
9
12
16

方式2：修正手动实现的模板函数

针对输入向量直接处理，自动计算输出维度，避免手动指定错误：

#include "Eigen/Dense"
#include <iostream>

using Index = Eigen::Index;

template <typename T, Index N>
constexpr inline Eigen::Vector<T, N*(N+1)/2>
get_upper_triangular_vec(const Eigen::Vector<T, N>& vec) {
    const auto outer_prod = vec * vec.transpose();
    Eigen::Vector<T, N*(N+1)/2> result;
    Index pos = 0;
    for (Index i = 0; i < N; ++i) {
        const Index segment_len = N - i;
        // 提取第i行从第i列开始的连续元素
        result.segment(pos, segment_len) = outer_prod.row(i).segment(i, segment_len);
        pos += segment_len;
    }
    return result;
}

int main() {
    Eigen::Vector<double, 4> o_lin;
    o_lin << 1, 2, 3, 4;

    auto o_nonlin = get_upper_triangular_vec(o_lin);
    std::cout << o_nonlin << std::endl;
}

关键注意事项 #

• Eigen的视图类（如triangularView）保留原矩阵维度，不能直接转为向量，需通过内存映射或手动提取元素。
• 对称矩阵的上三角提取优先用selfadjointView::packedView，内存连续且效率更高。
• 编译时固定大小的向量需确保维度计算正确（n维向量对应上三角元素数为n*(n+1)/2），不确定维度时可改用Eigen::Dynamic类型。

内容的提问来源于stack exchange，提问作者Ben Andrews