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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Tue, Apr 6, 2010 at 3:03 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> 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.

 
oops AnyMatrixBase has been renamed EigenBase ...

gael

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

Finaly, products and solving would idealy be implemented using the ReturnByValue stuff.

The same design is also used for quaternions. (Eigen/src/Geometry/Quaternion.h)

Hope that helps.

gael.
 

  For a Toeplitz matrix, it is necessary to store only the first row and
first column. Therefore, many arithmetic ops become simple.

Thanks,
Manoj