symbolic calculation - speeding up

[Arda]

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,

[Christopher Creutzig]

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