Improve running time of code with integrals

Hello,

I'm aiming to create the following covariance matrix,

,

where the set of possible combinations is given by

.

This results in different combinations, meaning Σ is of dimension .

The integrand is given by with some parameters .

Below is my created function, which takes as input an array of parameters , the output is the resulting covariance matrix Σ.

Because the covariance matrix is symmetric, I just compute the lower triangular part and add the upper triangular part in the last step.

Computing this takes very long, after 5 minutes I just got the first 200 lines of the matrix. Regarding that all this is just one step in computing a Likelihood-function, which I aim to maximize for numerically, the running time has to be reduced. The problem is obviously the part where I'm integrating in the middle of the code, can anyone give me an advice how to optimize this part?

function Sigma = covariance_matrix(P)

sig = @(t,x) (P(1) + P(2).*exp(P(3).*(x+t))) .* (P(4)+t) .* exp(-P(5).*t);

i = 0;
for x1 = 30:95              
    for t1 = 0:95-x1        
        i = i+1;
        j = 1;
        for x2 = 30:x1      
            t2 = 0;
            while j <= i && t2 <= 95-x2     
                Sigma(i,j) = integral(@(s) integral(@(u) sig(u,x1-1+s), ...
                    t1+1-s,t1+2-s).*integral(@(u) sig(u,x2-1+s), ...
                    t2+1-s,t2+2-s),0,1,'ArrayValued',true);
                t2 = t2+1;
                j = j+1;
            end
        end
    end
end
Sigma = Sigma + (Sigma-diag(diag(Sigma)))';

end

Answer 1

[Harald]

Hi,

it may be an alternative to calculate Sigma analytically using Symbolic Math Toolbox (please excuse any typos :) )

syms k c r a b tau x s t1 t2 u t x1 x2
sig(tau, s) = (k + c*exp(r*(tau+x)))*(a+t)*exp(-b*tau);
S = int(   int(sig(u, x1-1+s), u, t1+1-s, t1+2-s)  ... 
    .*  int(sig(u, x2-1+s), u, t2+1-s, t2+2-s), s, 0, 1)

Using matlabFunction, you can convert this to a function handle.

This is quite a long expression and may take a bit to evaluate as well, but an advantage I see is that you can likely evaluate this for several t2 in a vectorized operation.

Best wishes,

Harald

[MGee]

Many thanks, this helps a lot. The whole matrix is computed now in about three minutes. Below is my updated code if anyone is interested. If someone sees other things to optimize please let me know :)

function Sigma = covariance_matrix_syms(P)

syms k c r a b tau x s t1 t2 u x1 x2
sig(tau,x) = (k+c*exp(r*(tau+x)))*(a+tau)*exp(-b*tau)
sig1(u,s) = sig(u,x1-1+s)
sig2(u,s) = sig(u,x2-1+s)
int1 = int(sig1,u,t1+1-s,t1+2-s)
int2 = int(sig2,u,t2+1-s,t2+2-s)
S = int(int1*int2,s,0,1)
Sfh = matlabFunction(S)

i = 0;
for x1 = 30:95              
    for t1 = 0:95-x1        
        i = i+1;
        j = 1;
        for x2 = 30:x1      
            t2 = 0;
            while j <= i && t2 <= 95-x2     
                Sigma(i,j) = Sfh(P(1),P(2),P(3),P(4),P(5),t1,t2,x1,x2);
                t2 = t2+1;
                j = j+1;
            end
        end
    end
end
Sigma = Sigma + (Sigma-diag(diag(Sigma)))';

end

[Torsten]

The computation of Sfh should be done once and for all and not each time your function "covariance_matrix_syms" is called. Then use Sfh as a second input to the function.

[MGee]

Thank you Torsten, this is also a good advice to me. I will change that.

[MGee]

As Harald mentioned in his answer, using vectorized operation reduces computing time again. The covariance matrix and mean vector is now computed in one function, running time is 20 seconds.

syms k c r a b tau x s t1 t2 u x1 x2
sig(tau,x) = (k+c*exp(r*(tau+x)))*(a+tau)*exp(-b*tau);
sig1(u,s) = sig(u,x1-1+s);
sig2(u,s) = sig(u,x2-1+s);
int1 = int(sig1,u,t1+1-s,t1+2-s);
int2 = int(sig2,u,t2+1-s,t2+2-s);
int3 = int(sig1,u,0,t1+1-s);
Sig = int(int1*int2,s,0,1);
Mean = int(int1*int3,s,0,1);
Sfh = matlabFunction(Sig);
Mfh = matlabFunction(Mean);

function [Sigma,M] = covariance_mean(P,Mfh,Sfh)

M = zeros(2211,1);
Sigma = zeros(2211);
Comb = zeros(2211,2211,4);
i = 0;
for x1 = 30:95              
    for t1 = 0:95-x1        
        i = i+1;
        j = 1;
        for x2 = 30:x1      
            t2 = 0;
            while j <= i && t2 <= 95-x2     
                Comb(i,j,1) = t1;
                Comb(i,j,2) = t2;
                Comb(i,j,3) = x1;
                Comb(i,j,4) = x2;
                t2 = t2+1;
                j = j+1;
            end
        end
    end
end
Comb = Comb + tril(NaN(2211),-1)';
M = Mfh(P(1),P(2),P(3),P(4),P(5),Comb(:,1,1),Comb(:,1,3));
Sigma = Sfh(P(1),P(2),P(3),P(4),P(5),Comb(:,:,1),Comb(:,:,2), ...
        Comb(:,:,3),Comb(:,:,4));
Sigma(isnan(Sigma)) = 0;
Sigma = Sigma + tril(Sigma,-1)';
M = M + 0.5.*diag(Sigma);

end

[Dyuman Joshi]

Another modification can be to combine this

Comb(i,j,1) = t1;
Comb(i,j,2) = t2;
Comb(i,j,3) = x1;
Comb(i,j,4) = x2;

to

Comb(i,j,1:4) = [t1 t2 x1 x2];