title: Margin of optimisation in C++ armadillo code
body: <p>I am trying to migrate a quite complex matlab funciton in c++ using the library <code>armadillo</code> but I have serious problems in terms of performance (my c++ version is much slower than the matlab one which is a bit odd). I was wondering whether any of you can spot a point where I can improve my code and possibly can give me some suggestions. The problem is created by a bit in which I try to minimise a function. This is the matlab code</p>
  
  <pre><code>lambda = fminbnd(@SSE, lower_limit,upper_limit,options, y0, x); 
  
  function L=SSE(lambda, alpha, y) 
  N=size(y,2);
  len=size(y,1);
  for i=1:N
      z(:,i)=jglog(y(:,i),alpha,lambda);
  end
  s = 0;    
  mean_spec=mean(z,2);
  for i=1:N        
      for j=1:len            
          s = s + (z(j,i)-mean_spec(j,1))^2;
      end
  end    
  L=s;
  
  end
  
  function z=glog(y,alpha,lambda) % Glog transform  
  z=log((y-alpha)+sqrt((y-alpha).^2+lambda));
  end
  
  function [zj, gmn]=jglog(y,y0,lambda)    
  z=glog(y,y0,lambda);
  gmn=exp(mean(log(sqrt((y-y0).^2+lambda))));
  zj=z.*gmn;
  end
  </code></pre>
  
  <p>Following <a href="https://stackoverflow.com/questions/19751636/is-there-a-faster-version-of-fminbnd-in-matlab">this post</a> I downloaded the c++ version of the matlab minimisation code (brent's method). You can find it <a href="http://people.sc.fsu.edu/~jburkardt/cpp_src/brent/brent.cpp" rel="nofollow noreferrer" title="here">here</a>.</p>
  
  <p>This is my c++ version. </p>
  
  <pre><code>class SSE_c : public brent::func_base   //create functor
  {
      public:
          mat A;
          double offset;
          virtual double operator() (double lam)
          {
              return SSE(A,offset,lam);
          }
          SSE_c(mat a,double of) {A=a;offset=of;}
  };
  
  SSE_c fun(x,y0);
  
  brent::glomin(low_limit,up_limit,c,100,step_threshold,step_threshold,fun,lambda);
  
  double SSE(mat&amp; m,double ofs,double lam)
  {
      mat z(m.n_rows,m.n_cols);
  
      for(uint i=0;i&lt;m.n_cols;i++)
          z.col(i) = jglog(m.col(i),ofs,lam);
  
      std::cout &lt;&lt; "Iteration:" &lt;&lt; count++;
  
      double s=0;
      vec mean_spec(z.n_rows);
      FuncOnMatRows(z,mean_spec,[](rowvec const&amp; w){return mean(w);});
      for(uint i=0;i&lt;z.n_cols;i++)
          for(uint j=0;j&lt;z.n_rows;j++)
              if(is_finite(z(j,i)))
                  s += pow((z(j,i)-mean_spec(j)),2);
      return s;
  }
  
  vec jglog(vec&amp;&amp; v,double ofs,double lam)
  {
      vec g=glogF(v,ofs,lam);
      double gmn;
      vec interm = log(sqrt(square(v-ofs)+lam));
      if(interm.is_finite())
          gmn=exp(mean(interm));
      else
          gmn=exp(mean(interm.elem(find_finite(interm))));    
      g = g*gmn;    
      return g;
  }
  
  vec glogF(vec&amp;&amp; v,double ofs,double lam)
  {
      return glogF(v,ofs,lam);
  }
  
  vec glogF(vec&amp; v,double ofs,double lam)
  {
      vec z = log((v-ofs)+sqrt(square(v-ofs)+lam));
      return z;
  }
  
  template&lt;typename Func&gt;
  void FuncOnMatRows(const mat&amp; M,vec&amp; v,Func const &amp; func)
  {
      for(uint i=0;i&lt;M.n_rows;i++)     // operation calculated on rows
      {
          if(M.row(i).is_finite())
              v(i) = func(M.row(i));
          else if(!any(M.row(i)&gt;0))
              v(i) = NAN;
          else
          {
              rowvec b=M.row(i);
              v(i) = func(b.elem(find_finite(b)).t());    //output of .elem is always colvec, so transpose
          }
      }
  }
  </code></pre>