And if it is of any use to scope Eigen malloc checks, I have a class here:
https://github.com/eacousineau/amber_de ... ig.hpp#L92
in addition to a few simple unittests:
https://github.com/eacousineau/amber_de ... config.cpp

Note that the unittests for enabling EIGEN_RUNTIME_NO_MALLOC in one source file and disabling it in EIGEN_RUNTIME_NO_MALLOC in other source files fail for builds with debug symbols due to inlining behavior, at least with g++-4.6.3 on Ubuntu 12.04.

It's not a good idea to explicitly compute the inverse of a matrix, unless it is close to be unitary and very small. In general this is numerically less accurate and much more expensive than calling solve on the factorized matrix: computing the inverse first requires to factorize the matrix and then call solve with a matrix on the right-hand-side.

eacousineau wrote:That is quite awesome indeed!

I had also struggled with computing an inverse and caching it without allocation while also allowing it to be dynamically sized, so I used a combination of pointers and evalTo().

Code: Select all

// For preallocating inverters
typedef Eigen::PartialPivLU<Eigen::MatrixXd> MatrixInverter;
typedef boost::shared_ptr<MatrixInverter> MatrixInverterPtr;
...
class Test
{
  MatrixXd A;
  MatrixXd Ainv;
  MatrixInverterPtr inverter;

  void resize(int n)
  {
     A.resize(n, n);
     Ainv.resizeLike(A);
     inverter.reset(new MatrixInverter(n));
  }
 
  void update(...)
  {
    ...
    inverter->compute(A).inverse().evalTo(Ainv);
    ...
  }
};


The `update()` method passes the `set_is_malloc_allowed(false)` check without the need to allocate on the heap.
However, I guess there could then be issues with heap-fragmentation.

And regarding numericaly efficiency, is it more efficient to use the solver (especially for conditioning) as opposed to caching the matrix inverse?

wow, i'll try your code as soon as possible! looks great.

if you want to solve Ax = b, i read everywhere that taking the inverse wan't good at all for 1 stability, 2 speed.

but maybe in some cases it will work.

there are so many types of matrices that it depends i guess.

Jeff

okay so now i have another challenge

one word : ROOTS



here is some code i put up that doesn't raise the "malloc error"

Code: Select all

 Eigen::internal::set_is_malloc_allowed(true);

     const int maxPoly = 21;
     int sizea = 7; // seven coeffs from a0 to a6

          double * coeffs_an = new double[sizea];

        coeffs_an[0] = 0.0001;
         ...


     Eigen::VectorXd hardCase_polynomial(maxPoly);

     Eigen::internal::set_is_malloc_allowed(false);

     Eigen::PolynomialSolver<double, maxPoly-1> solver;
    
    
     for (int i = 0; i < siza; i++)
     {
        hardCase_polynomial(i)= coeffs_an[i];
     }
      for (int i = siza; i < maxPoly; i++)
     {
        hardCase_polynomial(i)= 0;
     }

   solver.compute(hardCase_polynomial);
   
   const Eigen::PolynomialSolver<double, maxPoly-1>::RootsType & r= solver.roots();
   //r = solver.roots().transpose();


   cout << r<<endl;

but the thing is, it doesn't work.

What i thought was :
if a polynomial is of the form a0 + a1 *x + a2 * x*x + a3 * x*x*x ... + an * x^n, we could zeros some coeffs from "am" to "an", and still be able to compute the right roots.

of course it would maybe take a little bit more computations, but in some cases, it's manageable (with nowadays CPU power)

the only problem is that when i put all coeffs to 0 after index sizea, the solver reduces the size of the poly automatically and raises an error:
Assertion failed : Scalar(0) != poly[poly.size()-1]

even when i try to trick the solver by putting a very small value for "an"

so...is there a way to compute roots in real-time? with my code or something else?

thanks

Jeff.

The Polynomials module is in unsupported/ meaning that it did not received much polishing. Feel free to propose a patch fixing this issue. It should probably be enough to compute the actual degree in Companion.h and use the template Degree parameter as a maximal size in teh declaration of matrix/vector objects.

i got an issue once again.

this time it's about a "simple" matrix multiplication

the matrices are 512x512.

I've succeeded in using multiplication with fixed-size matrices, but with matrices this big (512x512) it cannot be done this way.

so i thought i would perform my own multiplication algorithm. Only problem it's soooo slooow (512x512x512 multiplications each time!)

Code: Select all

void RNN::Multiply(MatrixXd &out, MatrixXd &in1, MatrixXd &in2)
{

   for (int i = 0; i < in1.rows(); i++)
   {
      for (int j = 0; j < in2.cols(); j++)
      {
         (out)(i, j) = 0;
         for (int k = 0; k < in2.rows(); k++)
         {
            (out)(i, j) += (in2)(k, j) * (in1)(i, k);
         }
      }
   }
}

So what are my options here ? Maybe separating the big matrices into smaller matrices ? but i don't know how to do that...

Are there other ways i wouldn't've thought of?

thanks for help.

Jeff

You should not have issue with our product implementation, however, do not forget the noalias() to avoid the creation of a temporary:

C.noalias() = A * B;

okay so finally i found how to perform block multiplications, it's super simple:

you just name some sub-matrices and do the mul like if it was scalars.

for example : A = [1 2 3 4; 3 4 5 6; 4 5 6 7; 5 6 7 8]
you put P1= [1 2; 3 4] p2 = [3 4 ; 56] ...etc
and put the mul as : A = [P1 P2; P3 P4] and do the mul normally.

i ended up with some code that works for Matrices with size >=N. You can change the value of N of course.
The higher N, the faster the computation

I did a quick benchmark.
in debug mode :
naïve computing : 48 sec
optimized N=128 : 12 sec
optimized N = 32 : 16 sec

in release mode, naïve computing is less than 1 sec, so it globally goes as fast as in Matlab (which was the goal)

ooops almost forgot the code:

Code: Select all

void Multiply_fast(MatrixXd &out, MatrixXd &in1, MatrixXd &in2)
{
   /*here we separate a big matrices into smaller ones.
   //The multiplication is the same. if we replace A and B by:
          P1 P2              R1 R2
   A2 = [ P3 P4 ] and B2 = [ R3 R4 ]

                  a1 a2   
         where P1 = [a3 a4] etc.

         */
   static const int N = 128;

   int na1 = in1.rows()/ N;
   int na2 = in1.cols()/ N;
   int nb1 = in2.rows()/ N;
   int nb2 = in2.cols()/ N;

   static Matrix<double, Dynamic, Dynamic, 0, N, N> multi_temp1;
   static Matrix<double, Dynamic, Dynamic, 0, N, N> multi_temp2;
   static Matrix<double, Dynamic, Dynamic, 0, N, N> multi_temp3;

   if (multi_temp1.size() <1)
   {
      multi_temp1.resize(N,N);
      multi_temp2.resize(N,N);
      multi_temp3.resize(N,N);
   }



   for (int i = 0; i < na1; i++)
   {
      for (int j = 0; j < nb2; j++)
      {
         for (int k = 0; k < nb1; k++)
         {
            for (int o = 0; o < N; o++)
            for (int p = 0; p < N; p++)
            {
               (multi_temp1)(o,p) = in1(i*N + o, j*N + p);
               (multi_temp2)(o,p) = in2(i*N + o, j*N + p);
            }

            multi_temp3 = multi_temp1 * multi_temp2;
            for (int o = 0; o < N; o++)
            for (int p = 0; p < N; p++)
            {
               (out)(i*N+o, j*N+p) = multi_temp3(o,p);
            }
         }
      }
   }
}

tell me what you think of this code, can it be optimized?

Jeff

EDI1: oh just saw the above post. will try noalias() and see the speed

Thanks bro you rock!.

Jeff

ggael wrote:You should not have issue with our product implementation, however, do not forget the noalias() to avoid the creation of a temporary:

C.noalias() = A * B;

eh, noalias() raises a no malloc flag!

what's going on?

Does your platform support alloca? In doubt try to compile with -DEIGEN_ALLOCA=alloca