Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
symbolic calculation - speeding up

Subject: symbolic calculation - speeding up

From: Arda

Date: 15 Nov, 2012 21:24:21

Message: 1 of 2

I am using Symbolic Toolbox, and facing serious speed issues.
Before describing the problem I am facing - here is my intention:

I want to create a matrix A (size n by n) whose i , j entries are defined as follows:

A(i,j) = (i out of j)* p^i * p ^(j-i) * (-1)^(j-i)

This A matrix has a symbolic variable p, but as the size grows the computation(high orders of p, and combination(j,i) ) becomes too slow > over 2 hours for a matrix of size 400*400.

And I would like to get a vector x = A*b
where b is an 1*n vector whose values are to be read from data. I defined b as a symbolic variable so that for different data files, I would simply plug in and compute A*b in terms of p.

Here is the code I am using:

%%%%%%%%%%%%%%%%
    %%%fill matrix
 
   n =400;
     %n = length(rhs);
    A = sym(zeros(n,n));
    sym p;

    for i=1:n,
        for j=i:n,

            a = j -i;
             
                
                if i <50
                            %Exact Binomial Probability
                            comb = exp( gammaln(i+a+1) - gammaln(i+1) - gammaln(a+1) );
                            A(i,j) = ((-1)^(a))*comb*((1-p)^(a));

                else
                            %%Poisson approximation
                            A(i,j) = ((-1)^(a))*exp(-(i+a)*p + i*log((i+a)*p) - 0.5*log(2*pi*i) + (i+a)*log(i/(exp(1))) )/(p^i);

    end
        end
    end
    
    
    b = sym('B',[1 n])
    b = sym(b, 'real')
    

         
   xvect = A * b';
    %%And at the end I want to able to plug in rhs and get the x vector in
    %%terms of p. After I will evaluate different p values to see which one
    %%does make sense in my problem.
    
   subs(xvect,sym('B',[1 n]),rhs);
   
  
   %%I haven't included all the part of binomial formula when computing A
   %%(to perfom less calculations)
   %%thus here I multiply them by these coefficients.
    post =sym(zeros(1,n));
 
    for i=1:n
       
        post(i)=p ^i;
        
    end
   xvect = xvect.*post';

%%%%%%%%%%%%%%%%


In a thread I'have read that using MatlabFunction is a better idea, but still the following operation : Z = matlabFunction(x_rhs,'vars',[p, sym('B',[1 n])])
takes hours to perfom.

I would appreciate any suggestions to speed up my code.

Best regards,

Subject: symbolic calculation - speeding up

From: Christopher Creutzig

Date: 3 Dec, 2012 09:29:54

Message: 2 of 2

On 15.11.12 22:24, Arda wrote:
> I am using Symbolic Toolbox, and facing serious speed issues.
> Before describing the problem I am facing - here is my intention:
>
> I want to create a matrix A (size n by n) whose i , j entries are defined as follows:
>
> A(i,j) = (i out of j)* p^i * p ^(j-i) * (-1)^(j-i)

I assume you mean (1-p)^(j-i)? Otherwise, the p factors combine to p^j.

> This A matrix has a symbolic variable p, but as the size grows the computation(high orders of p, and combination(j,i) ) becomes too slow > over 2 hours for a matrix of size 400*400.
>
> And I would like to get a vector x = A*b
> where b is an 1*n vector whose values are to be read from data. I defined b as a symbolic variable so that for different data files, I would simply plug in and compute A*b in terms of p.

That sounds numerically unstable. Have you checked that writing the
formula this way is actually a good idea on your data?

> Here is the code I am using:
>
> %%%%%%%%%%%%%%%%
> %%%fill matrix
>
> n =400;
> %n = length(rhs);
> A = sym(zeros(n,n));
> sym p;

That should be

p = sym('p');

since the sym command does not assign to a variable (and using syms
inside procedures may confuse the linter and possibly the optimizer).

> In a thread I'have read that using MatlabFunction is a better idea, but still the following operation : Z = matlabFunction(x_rhs,'vars',[p, sym('B',[1 n])])
> takes hours to perfom.

It may be worth looking at one of the entries you created. If it
contains very complicated constants, you may be better of switching to
vpa arithmetic early (i.e., when putting the constants into the
expressions in the first place), either by explicitly calling the vpa
function or by converting your doubles with one of the additional
arguments to sym, such as 'd'. If, on the other hand, your formulas look
more complicated than they need to be, add a call to simplify, probably
best after computing A*b'. Although with your particular formulas, I
don't expect that would do much.


HTH,
Christopher

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us