c++recurrent-neural-networkeigen3

Mapping Eigen triangularView to an Eigen vector


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, ononlin, from the linear feature vector, olin. For example to construct a quadratic nonlinear feature vector, the following procedure is followed:

In the paper, they give this operation the symbol: op

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.

Solution

  • 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