Mapping Eigen triangularView to an Eigen vector

Question

I'm currently attempting to implement a new version of a reservoir neural network (like a recurrent neural network but only the output layer weights are trained) according to it's description in this paper.

Part of this procedure describes how to construct a nonlinear feature vector, , from the linear feature vector, . For example to construct a quadratic nonlinear feature vector, the following procedure is followed:

• All monomials (singular term polynomial) of the quadratic polynomial are captured by the outer product
• The nonlinear feature vector is composed of the unique monomials which comprise the upper triangular elements of the outer product tensor

In the paper, they give this operation the symbol:

As an example, the following linear feature vector:

1
2
3
4

After calculating the outer product with itself becomes,

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

The nonlinear feature vector is then, by my understanding,

1
2
3
4
4
6
8
9
12
16

Ideally I'm looking for is a procedure using Eigen which allows me to do this?

For examples, I've tried something like this:

// You'll need to have the Eigen library
#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() );

    // This fails @ runtime for me
    upper_tri.evalTo( upper_tri_vec );
    std::cout << upper_tri_vec << std::endl;

    // This won't compile as the triangularView doesn't have a .reshaped() method
    // Eigen::Vector<double, 10> o_nonlin{ upper_tri.reshaped() };
}

Which gives the output:

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

Is there some method I can use to remedy this?

--- EDIT ---

While I've been messing around to try & implement this specific upper-triangular elements -> vector mapping operation I've run into a different issue when trying something like this:

// Shorter name for 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;
}

However, this gives the following runtime error:

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

Answer 1

I would use a Map (documentation here) to basically copy the triangularView to the Vector.

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

int main(int argc, char *argv[]) {

    Eigen::Matrix<double, 4, 1> olin;
    olin << 1,2,3,4;
    Eigen::Matrix4d outer_product;
    outer_product = olin * olin.transpose();
    Eigen::Vector<double, 16> triview_vec;
    Eigen::Matrix<double, 4, 4>::Map(triview_vec.data(), 16, 1) = outer_product.triangularView<Eigen::Upper>();
    std::cout << triview_vec << std::endl;
    return 0;
}

At this point triview_vec will look like:

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

Then you can use whatever method you like best to pull out the values you care about. I am not good at indexing into Eigen vectors so here is my brute force way of doing it.

Eigen::Vector<int, 16> down = (triview_vec.array() > 0.0).select(Eigen::Vector<int, 16>::Ones(),0);
Eigen::Vector<double, Eigen::Dynamic> nonlin;
int size_nonlin = down.sum();
nonlin.resize(size_nonlin);

size_t tracker = 0;
for (size_t i = 0; i < 16; ++i){
    if (down(i)){
        nonlin(tracker) = triview_vec(i);
        tracker++;
    }
}

The end result is:

 1
 2
 4
 3
 6
 9
 4
 8
12
16