Re: [eigen] Efficient syntax for Ax=b where A is diagonal/Toeplitz

[Manoj Rajagopalan <rmanoj@xxxxxxxxx>]

Hi Gael (and other developers),

   A quick question regarding a Toeplitz matrix implementation based on 
pointers you'd mentioned earlier ...


On Tuesday 06 April 2010 09:03:09 am Gael Guennebaud wrote:
> On Fri, Apr 2, 2010 at 9:32 PM, Manoj Rajagopalan <rmanoj@xxxxxxxxx> wrote:
> > On Friday 02 April 2010 02:21:59 pm Gael Guennebaud wrote:
> > > On Fri, Apr 2, 2010 at 5:03 PM, Manoj Rajagopalan <rmanoj@xxxxxxxxx>
> >
> > wrote:
> > > >    Also, efficient storage, arithmetic and linear-system solution
> > > > (O(n^2)) algorithms are possible when A is (symmetric) Toeplitz
> > > > (occur frequently in signal-processing and Fourier-series solution of
> >
> > ODE/PDE).
> >
> > > > Does the roadmap consider this case?
> > >
> > > I'd say contributions are welcome ;)
> > >
> > > gael
> >
> > Any hints on how to go about this? I could try my hand - I am experienced
> > C++
> > developer (yes, expr templates, CRTP/barton-nackman and all the related
> > good
> > stuff also). Just that I just got to know of Eigen. What to inherit? How
> > and
> > where to introduce specialized arithmetic functors into the existing
> > framework?
>
> I would implement it in the same vein as DiagonalMatrix/DiagonalWrapper:
> - add a ToeplitzBase class inheriting AnyMatrixBase. In addition to the
> rows() and cols() functions, all classes inheriting ToeplitzBase should
> provide a coeffs() function returning a CoeffType reference to the dense
> vector expression representing the first row and first column. ToeplitzBase
> will implement the common API for all Toeplitz objects.
> - add two classes ToeplitzMatrix & ToeplitzWrapper inheriting ToeplitzBase.
> The former would store a single dense vector + rows and cols with the
> appropriate functions allowing to manipulate the matrix. The later would
> simply store a reference to a dense expression representing the first
> column and first row. Then, e.g., the addition of two Toeplitz objects
> would be implemented in ToeplitzBase by nesting the expression of the
> addition of the two dense vectors into a ToeplitzWrapper:
>
> template<typename OtherDerived>
> ToeplitzWrapper<CwiseBinaryOp<ei_scalar_add_op<Scalar>, CoeffType, typename
> OtherDerived::CoeffType>
> operator+(const ToeplitzBase<OtherDerived>& other)
> {
>   return ToeplitzWrapper<CwiseBinaryOp<ei_scalar_add_op<Scalar>, CoeffType,
> typename OtherDerived::CoeffType>(this->coeffs() + rhs.coeffs());
> }
>

At the assign stage, when the CwiseBinaryOp<> expression is being referenced, 
it seems this object has not being constructed properly. My debugger tells me 
it has difficulty reading the CwiseBinaryOp<> member of ToeplitzWrapper<>. 
Could this be because it is a reference which was assigned to an instance 
created at the '+' but is no longer alive?

My implementation is based on DiagonalMatrix and DiagonalWrapper as you 
mentioned.

I am attaching my .cpp file and a test program in case you want to take a 
look. With DiagonalMatrix, no such addition operation is provided.

One alterative is to not allow such addition operations for Toeplitz in which 
case, one would have to resort to the roundabout method I have used in 
the .cpp test file.

BTW, once I complete developing this class, do you think it might be a useful 
addition to Eigen3, perhaps among the unsupported modules?

thanks,
Manoj

#include <Eigen/Core?
#include <SymmetricToeplitzMatrix.h>
#include <iostream>

using namespace Eigen;
using namespace std;

int main(void)
{
	int const N = 5;
	SymmetricToeplitzMatrixXi S1(N), S2(N), S3(N);
	S1.coeffVector().setRandom();
	S2 = S1;
	cout << "S2 =\n" << S2.toDenseMatrix() << endl;
    
    // the following does NOT work
    // S3 = S1 + S2;
    
    // the following works instead
    S3 = SymmetricToeplitzWrapper<CwiseBinaryOp<ei_scalar_sum_op<int>,VectorXi,VectorXi> >(S1.coeffVector() + S2.coeffVector());
	
    cout << "S3 =\n" << S3.toDenseMatrix() << endl;
	return 0;
}

#ifndef EIGEN_SYMMETRIC_TOEPLITZ_MATRIX_H
#define EIGEN_SYMMETRIC_TOEPLITZ_MATRIX_H

namespace Eigen {

// Forward declarations
template<typename _Scalar, int _Rows, int _MaxRows>
class SymmetricToeplitzMatrix;

template<typename _CoeffVectorType>
class SymmetricToeplitzWrapper;



template<typename Derived>
class SymmetricToeplitzBase : public EigenBase<Derived>
{
public:
    typedef typename ei_traits<Derived>::CoeffVectorType CoeffVectorType;
    typedef typename CoeffVectorType::Scalar Scalar;
    typedef typename ei_traits<Derived>::StorageKind StorageKind;
    typedef typename ei_traits<Derived>::Index Index;
    enum {
        RowsAtCompileTime = CoeffVectorType::RowsAtCompileTime,
        ColsAtCompileTime = RowsAtCompileTime,
        MaxRowsAtCompileTime = CoeffVectorType::MaxRowsAtCompileTime,
        MaxColsAtCompileTime = MaxRowsAtCompileTime,
        IsVectorAtCompileTime = 0,
        Flags = 0
    };
    
    typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType;
    
    inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
    inline Derived& derived() { return *static_cast<Derived*>(this); }
    
    DenseMatrixType toDenseMatrix() const { return derived(); }
    template<typename DenseDerived>
        void evalTo(MatrixBase<DenseDerived> &other) const;
    template<typename DenseDerived>
        void addTo(MatrixBase<DenseDerived> &other) const;
    template<typename DenseDerived>
        void subTo(MatrixBase<DenseDerived> &other) const;
    
    inline const CoeffVectorType& coeffVector() const { return derived().coeffVector(); }
    inline CoeffVectorType& coeffVector() { return derived().coeffVector(); }
    
    inline Index rows() const { return coeffVector().size(); }
    inline Index cols() const { return coeffVector().size(); }
    
    /** Addition */
    template<typename OtherDerived>
        SymmetricToeplitzBase<SymmetricToeplitzWrapper<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,
                                                                     CoeffVectorType,
                                                                     typename OtherDerived::CoeffVectorType>
                                                      >
                             > operator+(const SymmetricToeplitzBase<OtherDerived>& other)
    { return SymmetricToeplitzWrapper<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,
                                                    CoeffVectorType,
                                                    typename OtherDerived::CoeffVectorType>
                                     >(this->derived().coeffVector() + other.derived().coeffVector());
    }
    
    /** Subtraction */
    template<typename OtherDerived>
        SymmetricToeplitzBase<SymmetricToeplitzWrapper<CwiseBinaryOp<ei_scalar_difference_op<Scalar>,
                                                                     CoeffVectorType,
                                                                     typename OtherDerived::CoeffVectorType>
                                                      >
                             > operator-(const SymmetricToeplitzBase<OtherDerived>& other)
    { return SymmetricToeplitzWrapper<CwiseBinaryOp<ei_scalar_difference_op<Scalar>,
                                                    CoeffVectorType,
                                                    typename OtherDerived::CoeffVectorType>
                                     >(this->coeffVector() - other.coeffVector());
    }

    /** @todo ToeplitzProduct */
//     template<typename MatrixDerived>
//         const ToeplitzProduct<MatrixDerived, Derived, OnTheLeft>
//         operator*(const MatrixBase<MatrixDerived> &matrix) const;
    
    /** @todo inverse (also a SymmetricToeplitzMatrix) */
};


template<typename Derived>
template<typename DenseDerived>
void SymmetricToeplitzBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
{
    other.diagonal().setConstant(coeffVector()[0]);
    for(Index diag = 1; diag < rows(); ++diag) {
        other.diagonal(diag).setConstant(coeffVector()[diag]);
        other.diagonal(-diag).setConstant(coeffVector()[diag]);
    }
}


template<typename Derived>
template<typename DenseDerived>
void SymmetricToeplitzBase<Derived>::addTo(MatrixBase<DenseDerived> &other) const
{
    other.diagonal().array() += coeffVector()[0];
    for(Index diag = 1; diag < rows(); ++diag) {
        other.diagonal(diag).array() += coeffVector()[diag];
        other.diagonal(-diag).array() += coeffVector()[diag];
    }
}


template<typename Derived>
template<typename DenseDerived>
void SymmetricToeplitzBase<Derived>::subTo(MatrixBase<DenseDerived> &other) const
{
    other.diagonal().array() -= coeffVector()[0];
    for(Index diag = 1; diag < rows(); ++diag) {
        other.diagonal(diag).array() -= coeffVector()[diag];
        other.diagonal(-diag).array() -= coeffVector()[diag];
    }
}



template<typename _CoeffVectorType>
struct ei_traits<SymmetricToeplitzWrapper<_CoeffVectorType> >
{
    typedef _CoeffVectorType CoeffVectorType;
    typedef typename CoeffVectorType::Scalar Scalar;
    typedef typename CoeffVectorType::Index Index;
    typedef typename CoeffVectorType::StorageKind StorageKind;
    enum {
        RowsAtCompileTime = CoeffVectorType::RowsAtCompileTime,
        ColsAtCompileTime = RowsAtCompileTime,
        MaxRowsAtCompileTime = CoeffVectorType::MaxRowsAtCompileTime,
        MaxColsAtCompileTime = MaxRowsAtCompileTime,
        Flags = 0
    };
};


template<typename _CoeffVectorType>
class SymmetricToeplitzWrapper
: public SymmetricToeplitzBase<SymmetricToeplitzWrapper<_CoeffVectorType> >,
    ei_no_assignment_operator
{
    typedef SymmetricToeplitzBase<SymmetricToeplitzWrapper> Base;
public:
    #ifndef EIGEN_PARSED_BY_DOXYGEN
    typedef _CoeffVectorType CoeffVectorType;
    typedef SymmetricToeplitzWrapper Nested;
    typedef typename Base::DenseMatrixType::Index Index;
    #endif
    
    /** Constructor from expression of coefficients to wrap. */
    inline SymmetricToeplitzWrapper(const CoeffVectorType& coeffVec) : m_coeffVector(coeffVec) {}
    
    /** \returns a const reference to the wrapped expression of coefficients. */
    const CoeffVectorType& coeffVector() const { return m_coeffVector; }
    
    inline Index rows() const { return m_coeffVector.size(); }
    inline Index cols() const { return m_coeffVector.size(); }
    
protected:
    const typename CoeffVectorType::Nested m_coeffVector;
};


template<typename _Scalar, int _Rows, int _MaxRows>
    struct ei_traits<SymmetricToeplitzMatrix<_Scalar,_Rows,_MaxRows> >
{
    typedef _Scalar Scalar;
    typedef Matrix<_Scalar,_Rows,1,0,_MaxRows,1> CoeffVectorType;
    typedef typename CoeffVectorType::StorageKind StorageKind;
    typedef typename CoeffVectorType::Index Index;
};


template<typename _Scalar, int _Rows, int _MaxRows>
class SymmetricToeplitzMatrix
: public SymmetricToeplitzBase<SymmetricToeplitzMatrix<_Scalar,_Rows,_MaxRows> >
{
public:
    #ifndef EIGEN_PARSED_BY_DOXYGEN
    typedef typename ei_traits<SymmetricToeplitzMatrix>::CoeffVectorType CoeffVectorType;
    typedef const SymmetricToeplitzMatrix& Nested;
    typedef _Scalar Scalar;
    typedef typename ei_traits<SymmetricToeplitzMatrix>::StorageKind StorageKind;
    typedef typename ei_traits<SymmetricToeplitzMatrix>::Index Index;
    #endif
    
protected:
    CoeffVectorType m_coeffVector;
    
public:
    /** const version of coeffVector(). */
    inline const CoeffVectorType& coeffVector() const { return m_coeffVector; }
    /** \returns a reference to the stored vector of coefficients. */
    CoeffVectorType& coeffVector() { return m_coeffVector; }
    
    /** Default constructor without initialization */
    inline SymmetricToeplitzMatrix() {}
    
    /** Constructs a diagonal matrix with given dimension  */
    inline SymmetricToeplitzMatrix(Index dim) : m_coeffVector(dim) {}
    
    /** 2D constructor. */
    inline SymmetricToeplitzMatrix(const Scalar& x, const Scalar& y) : m_coeffVector(x,y) {}
    
    /** 3D constructor. */
    inline SymmetricToeplitzMatrix(const Scalar& x, const Scalar& y, const Scalar& z) : m_coeffVector(x,y,z) {}
    
    /** Copy constructor. */
    template<typename OtherDerived>
        inline SymmetricToeplitzMatrix(const SymmetricToeplitzBase<OtherDerived>& other) : m_coeffVector(other.coeffVector()) {}
    
    #ifndef EIGEN_PARSED_BY_DOXYGEN
    /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */
    inline SymmetricToeplitzMatrix(const SymmetricToeplitzMatrix& other) : m_coeffVector(other.coeffVector()) {}
    #endif
    
        /** generic constructor from expression of the diagonal coefficients */
    template<typename OtherDerived>
        explicit inline SymmetricToeplitzMatrix(const MatrixBase<OtherDerived>& other) : m_coeffVector(other)
    {}
    
        /** Copy operator. */
    template<typename OtherDerived>
        SymmetricToeplitzMatrix& operator=(const SymmetricToeplitzBase<OtherDerived>& other)
    {
        m_coeffVector = other.coeffVector();
        return *this;
    }
    
    #ifndef EIGEN_PARSED_BY_DOXYGEN
    /** This is a special case of the templated operator=. Its purpose is to
      * prevent a default operator= from hiding the templated operator=.
      */
    SymmetricToeplitzMatrix& operator=(const SymmetricToeplitzMatrix& other)
    {
        m_coeffVector = other.m_coeffVector;
        return *this;
    }
    #endif
    
    /** Resizes to given size. */
    inline void resize(Index size) { m_coeffVector.resize(size); }
    /** Sets all coefficients to zero. */
    inline void setZero() { m_coeffVector.setZero(); }
    /** Resizes and sets all coefficients to zero. */
    inline void setZero(Index size) { m_coeffVector.setZero(size); }
    /** Sets this matrix to be the identity matrix of the current size. */
    inline void setIdentity() { m_coeffVector.setZero(); m_coeffVector[0] = Scalar(1); }
    /** Sets this matrix to be the identity matrix of the given size. */
    inline void setIdentity(Index size) { m_coeffVector.setZero(size); m_coeffVector[0] = Scalar(1); }
};


typedef SymmetricToeplitzMatrix<int,1,1> SymmetricToeplitzMatrix1i;
typedef SymmetricToeplitzMatrix<float,1,1> SymmetricToeplitzMatrix1f;
typedef SymmetricToeplitzMatrix<double,1,1> SymmetricToeplitzMatrix1d;
typedef SymmetricToeplitzMatrix<std::complex<float>,1,1> SymmetricToeplitzMatrix1cf;
typedef SymmetricToeplitzMatrix<std::complex<double>,1,1> SymmetricToeplitzMatrix1cd;

typedef SymmetricToeplitzMatrix<int,2,2> SymmetricToeplitzMatrix2i;
typedef SymmetricToeplitzMatrix<float,2,2> SymmetricToeplitzMatrix2f;
typedef SymmetricToeplitzMatrix<double,2,2> SymmetricToeplitzMatrix2d;
typedef SymmetricToeplitzMatrix<std::complex<float>,2,2> SymmetricToeplitzMatrix2cf;
typedef SymmetricToeplitzMatrix<std::complex<double>,2,2> SymmetricToeplitzMatrix2cd;

typedef SymmetricToeplitzMatrix<int,3,3> SymmetricToeplitzMatrix3i;
typedef SymmetricToeplitzMatrix<float,3,3> SymmetricToeplitzMatrix3f;
typedef SymmetricToeplitzMatrix<double,3,3> SymmetricToeplitzMatrix3d;
typedef SymmetricToeplitzMatrix<std::complex<float>,3,3> SymmetricToeplitzMatrix3cf;
typedef SymmetricToeplitzMatrix<std::complex<double>,3,3> SymmetricToeplitzMatrix3cd;

typedef SymmetricToeplitzMatrix<int,4,4> SymmetricToeplitzMatrix4i;
typedef SymmetricToeplitzMatrix<float,4,4> SymmetricToeplitzMatrix4f;
typedef SymmetricToeplitzMatrix<double,4,4> SymmetricToeplitzMatrix4d;
typedef SymmetricToeplitzMatrix<std::complex<float>,4,4> SymmetricToeplitzMatrix4cf;
typedef SymmetricToeplitzMatrix<std::complex<double>,4,4> SymmetricToeplitzMatrix4cd;

typedef SymmetricToeplitzMatrix<int,5,5> SymmetricToeplitzMatrix5i;
typedef SymmetricToeplitzMatrix<float,5,5> SymmetricToeplitzMatrix5f;
typedef SymmetricToeplitzMatrix<double,5,5> SymmetricToeplitzMatrix5d;
typedef SymmetricToeplitzMatrix<std::complex<float>,5,5> SymmetricToeplitzMatrix5cf;
typedef SymmetricToeplitzMatrix<std::complex<double>,5,5> SymmetricToeplitzMatrix5cd;

typedef SymmetricToeplitzMatrix<int,Dynamic,Dynamic> SymmetricToeplitzMatrixXi;
typedef SymmetricToeplitzMatrix<float,Dynamic,Dynamic> SymmetricToeplitzMatrixXf;
typedef SymmetricToeplitzMatrix<double,Dynamic,Dynamic> SymmetricToeplitzMatrixXd;
typedef SymmetricToeplitzMatrix<std::complex<float>,Dynamic,Dynamic> SymmetricToeplitzMatrixXcf;
typedef SymmetricToeplitzMatrix<std::complex<double>,Dynamic,Dynamic> SymmetricToeplitzMatrixXcd;

} // namespace Eigen

#endif // EIGEN_SYMMETRIC_TOEPLITZ_MATRIX_H