How can I reduce the execution time of my code?

I wrote this code. it is working for m=32 and I am getting output with in 2 minutes. But if I am putting m=64, it is running. I checked for 24 hours. I am not getting the output. Can you tell me how I can reduce the execution time of this code.

clc;
clear all;
close all;
a=1/3;
delta=1.3;
Gamma=0.6;
global m t matrix matrix1
m=32;
matrix = zeros(m, m);
matrix1=zeros(m,m);
for i=1:m
    matrix(1,i)= 1/sqrt(m);
    for l=1:m
        t(l)=(l-0.5)./m;
    end
end
for i = 1:m-1
    for j = 0:(i - 1)
        for k = 1:(2^j + 1)-1
            if i == 2^j + k - 1
                fprintf('For i = %d, j = %d, and k = %d\n', i, j, k)
                for l = 1:m % Changed variable name from j to l
                    t_val = t(l);
                    matrix1(0+1,l) = (t_val^a)/(gamma(1+a)*sqrt(m));
                    if (k - 1) / 2^j <= t_val
                        if t_val < (k - 0.5) / 2^j
                            matrix(i+1, l) = 1/sqrt(m) *(2^(j/2));
                            matrix1(i+1, l) = 2^(j/2) *(1/(gamma(1+a)*sqrt(m)) * (t_val - (k - 1) / 2^j)^a);
                        elseif t_val < k / 2^j
                            matrix(i+1, l) = 1/sqrt(m) * (-2^(j/2));
                            matrix1(i+1, l) = 2^(j/2) * (1/(gamma(1+a)*sqrt(m)) * (t_val - (k - 1) / 2^j)^a - 2 / (gamma(1+a)*sqrt(m)) * (t_val - (k - 0.5) / 2^j)^(a));
                            
                        elseif t_val < 1
                            matrix1(i+1, l) = 2^(j/2) * (1/(gamma(1+a)*sqrt(m))* (t_val - (k - 1) / 2^j)^(a) - (2 / (gamma(1+a)*sqrt(m)) * (t_val - (k - 0.5) / 2^j)^(a)) + (1 / (gamma(1+a)*sqrt(m)) * (t_val - k / 2^j)^(a)))
                            
                        else
                            matrix(i+1, l) = 0;
                            matrix1(i+1, l) = 0;
                        end
                    end
                end
            end
        end
    end
end
disp(matrix) % Display the resulting matrix
disp(matrix1) % Operational matrix
cstart=zeros(1,2*m);
G=fsolve(@nddt,cstart);
G
for i=1:m
    G1(i)=G(i);
end
for i=1:m
    u(i)=delta+G1*matrix1(:,i);
end
plot(t,u,'-b')
hold on
p=1;
for j=m+1:2*m
    G2(p)=G(j);
    p=p+1;
end
for i=1:m
    v(i)=Gamma+G2*matrix1(:,i);
end
plot(t,v,'-r')




function [F]=nddt(E)
delta = 1.3;
Gamma = 0.6;
global m t matrix matrix1
for l=1:m
    F1(l)=E(1:m)* matrix(:,l) - t(l)*E(1:m)* matrix1(:,1)- t(l)*delta +(E(1:m)* matrix1(:,1) + delta)*(E(m+1:2*m)* matrix1(:,1) + Gamma);
    F2(l)=E(m+1:2*m)* matrix(:,l) - (E(1:m) * matrix1(:,1) + delta)*(E(m+1:2*m)*matrix1(:,1)+ Gamma) + t(l)*(E(m+1:2*m) * matrix1(:,1) + Gamma);
end
F=[F1,F2];
end

[Stephen23]

Look at how your function scales:

m = pow2(0:5)

m = 1x6

1 2 4 8 16 32

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

v = uint64(arrayfun(@myfun,m))

v = 1x6

0 2 44 1976 1048304 137438952416

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

semilogy(v) % ouch

Now estimate/calculate how many iterations your code would perform for m=64. Then think about how how long each iteration takes, and thus how long your function would need to run for, in order to perform all of those iterations.

Not all problems can be solved using the naive approach of iterating though all permutations.

function n = myfun(m)
n = 0;
for i = 1:m-1
    for j = 0:(i - 1)
        for k = 1:(2^j + 1)-1
            n = n+m;
        end
    end
end
end

[Torsten]

OP wrote that the generation of the matrices is not the problem, but the subsequent call to "fsolve".

[Surath Ghosh]

G=fsolve(@nddt,cstart); This is the problem. For this comments, output is not coming. For m=64, basically using fsolve, I want to solve 64*2=128 equations.

[Walter Roberson]

for i = 1:m-1
    for j = 0:(i - 1)
        for k = 1:(2^j + 1)-1
            if i == 2^j + k - 1

That can be rewritten,

for i = 1:m-1
    for j = 0:(i - 1)
        k = i - 2^j + 1
        if k >= 1 & k <= 2^j

[Surath Ghosh]

Thank you very much sir for your comments. But I am getting the output that are the matrix and matrix1 within 10 minutes. That 4 loops are used only for matrix and matrix1. Main problem is in the G=fsolve(@nddt,cstart), when I called the nddt function. For this comment, output is not coming. For m=64, basically using fsolve, I want to solve 64*2=128 equations. The code is stucked. I tried to run this code for 96 hours but I did not get the output.

Answer 1

[Dinesh]

Hi Surath,

Here are some optimizations that I could identify in your code:

• Precompute Reusable Values: Avoid repetitive computation of constants or expressions within loops.

codegammaFactor = 1 / (gamma(1 + a) * sqrt(m));

• Vectorize Conditional Assignments: Use logical indexing to apply conditions across a vector or matrix.

codecondition = (k - 1) / 2^j <= t & t < k / 2^j;
matrix1(i, condition) = 2^(j/2) * gammaFactor .* (t(condition) - (k - 1) / 2^j).^a;

• Optimize Matrix Operations: For operations applied to each row or column, use matrix multiplication or dot products.

codeG1 = G(1:m);
u = delta + G1 * matrix1;

Further Optimization Strategies:

• Parallel Computing: For large-scale matrix operations or loops that can run independently, consider using MATLAB's parallel computing capabilities (parfor loops, distributed arrays).
• Profiling and Algorithm Refinement: Utilize MATLAB's profiler to identify bottlenecks. In some cases, rethinking the algorithm or applying mathematical simplifications can lead to significant performance gains.

Reference Links:

https://www.mathworks.com/help/parallel-computing/getting-started-with-parallel-computing-toolbox.html

https://www.mathworks.com/help/matlab/matlab_prog/profiling-for-improving-performance.html

[Surath Ghosh]

Thank you very much for your comments. I am getting matrix and matrix1 within 10 minutes. So, it is not taking more time. But using the below code, I am not getting the output (checked for 24 hours). When I called the nddt function, I am not getting the output.

function [F]=nddt(E)

delta = 1.3;

Gamma = 0.6;

global m t matrix matrix1

for l=1:m

F1(l)=E(1:m)* matrix(:,l) - t(l)*E(1:m)* matrix1(:,1)- t(l)*delta +(E(1:m)* matrix1(:,1) + delta)*(E(m+1:2*m)* matrix1(:,1) + Gamma);

F2(l)=E(m+1:2*m)* matrix(:,l) - (E(1:m) * matrix1(:,1) + delta)*(E(m+1:2*m)*matrix1(:,1)+ Gamma) + t(l)*(E(m+1:2*m) * matrix1(:,1) + Gamma);

end

F=[F1,F2];

end

[Torsten]

You could try "lsqnonlin" or "fmincon" instead of "fsolve" to solve your quadratic system of equations, but these are the only options you have.

Answer 2

[Walter Roberson]

for i = 1:m-1
    for j = 0:(i - 1)
        for k = 1:(2^j + 1)-1
            if i == 2^j + k - 1
                fprintf('For i = %d, j = %d, and k = %d\n', i, j, k)
                for l = 1:m % Changed variable name from j to l

You have 4 nested loops, three of which effectively have an upper limit of m, and the fourth has a much higher upper limit. When you double your m you have to expect that it will take notably longer.

[Surath Ghosh]

Thank you very much sir for your comments. But I am getting the output that are the matrix and matrix1 within 10 minutes. That 4 loops are used only for matrix and matrix1. Main problem is in the G=fsolve(@nddt,cstart), when I called the nddt function. For this comment, output is not coming. For m=64, basically using fsolve, I want to solve 64*2=128 equations. The code is stucked. I tried to run this code for 96 hours but I did not get the output.