Vectorization of for loops with symbolic functions inside

Hi,

I have been using for-loops extensively and they served the prupose well until I had to run a code that requires over 140 thousand iteration. At that point my computer started taking an eternity to process it. So, I tried searching for ways to optimize the code and something I got introduced to is vectorization of the for-loops.

I have tried looking at ways to vectorize the two for loops but I can't seem to figure out a way that would suit my case. In my code, inside each iteration I am using symbolic variables and solving them for their values. Have a look at the simplified version of my code below:

syms t
vec1 = [5*cosd(t); 5*sind(t)]; 
vec2 = [t^2+2; t^3+5];

start_time = 0;
end_time = 360;

for time1 = start_time:1:end_time-1
    for time2 = time1+1:1:end_time
        A = double(subs(vec1, t, time1));
        B = double(subs(vec2, t, time2));
        
        syms x y %two unknowns
        eq1 = transpose(A)*[x;y]
        eq2 = transpose(B)*[x;y]
        eq3 = [eq1; eq2];
        solutions = solve(eq3,[x,y]);
        
        %Output the solutions
        x.solution
        y.solution
    end 
end

Is there a way to vectorize the two time loops?

Thanks in advance for your help.

Answer 1

[Walter Roberson]

Yes, you can vectorize the code completely by just assigning 0 and 0 as the solutions.

You construct eq3 as the sum of two scalar terms each of which has x times something plus y times something. collect x and y to get something times x plus something times y. You solve the scalar sum for x and y (one equation with two unknown). But that has the trivial solution x=0 y=0 no matter what the details of the terms are.

[Mohammad Zainullah Khan]

My bad, I actually made a small mistake when writing it out here in the forum.

Eq3 is actually a 2by1 vector so it got two equations with two unknowns. Do you think I can somehow vectorize the two for loops of time1 and time2?

[Walter Roberson]

Yes, the solution is exactly the same.

eq1 is a1*x + a2*y

eq2 is b1*x + b2*y

by inspection, x=0 y=0 is a solution to the simultaneous equations.

[Mohammad Zainullah Khan]

Yeah, but that is for this particular code which is just a very simplified version that I posted to explain what I was looking for.

Let's say if the equations are a bit more challenging than that, for e.g

eq1 is a1*x + a2*y + a3

eq2 is b1*x + b2*y + b3

So, this would imply that there has to be only solution (or zero if lines are parallel). For this case, can I vectorize the two time for loops to improve time?

By vectorization, I mean something like this :link

[Torsten]

Please describe the problem that you really have, not a simplified version.

The two equations above can easily be solved analytically by Cramer's rule - and an analytical solution is the best vectorization.

[Mohammad Zainullah Khan]

My main query is related to this part

for time1 = start_time:1:end_time-1
    for time2 = time1+1:1:end_time
         
        %some formulas and equation solving here
        
    end 
end

In this part, instead of these nested for loops, can I can use vectorization such as done here (link). If yes, how? Sorry for any confusion earlier.

[Walter Roberson]

no you cannot.

Instead you solve the generalized problem in terms of t, and then either drop in the entire list of times, or loop (non-vectorized) using the precomputed formula.

It is a contradiction to talk about vectorization while requiring loops still exist.

[Torsten]

Maybe. But if you really have to solve some set of equations for all combinations of (ti,tj) for j>i, then your loop is the simplest way to do so. And usually, the loop will not be the time-consuming factor, but solving the equations therein.

[Mohammad Zainullah Khan]

Walter, the idea that you put forward of solving the generalized problem in terms of t and then dropping in the values of t at the end seems interesting. There may be cases where I can get away with it but generally for the complex problems that I am dealing with, it might mess up a lot of stuff.

[Mohammad Zainullah Khan]

Torsten, thank you for such a starightforward comment that the loop seems to be the simplest way to go.

Is there any optimizations that I can do when either using substitutions using 'subs' or while using the 'solve' function?

[Torsten]

For what final problem ?

As I said, solving

a1*x + a2*y + a3 = 0
b1*x + b2*y + b3 = 0

for x and y can be done analytically which will be way faster than "solve".

[Walter Roberson]

        A = double(subs(vec1, t, time1));
        B = double(subs(vec2, t, time2));

Instead before the loop

syms t1 t2
A = subs(vec1, t, t1) ;
B = subs(vec2, t, t2);

and proceed to create the equations and solve() for x and y. This will give you an analytical solution in terms of t1 t2

Now ndgrid to create T1 T2 from start_time:1:end_time. And subs those arrays into the x and y solutions for t1 and t2. You get out a full array for X and Y

You only want the triu() or tril() of those. It is not clear what you want to do with the solutions once you have them.

[Mohammad Zainullah Khan]

Torsten, what do you mean by solving it analytically? Do keep in mind that I would need to perform this for thousands of different equations and some might not have solutions (two lines parallel to each other).

[Mohammad Zainullah Khan]

Walter, I tried changing the code to the following as you suggested:

clear; close all; clc
syms t
vec1 = [5*cosd(t)+2; 5*sind(t)+1]; 
vec2 = [t^2+2; t^3+5];
start_time = 0;
end_time = 360;

syms t1 t2
A = subs(vec1, t, t1) ;
B = subs(vec2, t, t2);
syms x y %two unknowns
eq1 = transpose(A)*[x;y];
eq2 = transpose(B)*[x;y];
eq3 = [eq1+1; eq2+2];
solutions = solve(eq3,[x,y]);

[T1 T2] = ndgrid(start_time:60:end_time)
subs(solutions.x,[t1,t2],[T1 T2])
subs(solutions.y,[t1,t2],[T1 T2])

I am still getting an error on the last two line. I think I am not following the syntax correctly (quite possible since I haven't used ndgrid before).

[Walter Roberson]

Thousands of different equations? Or the same pattern of equations with thousands of different times?

If you use solve() to find the general pattern, then use numden() to split into numerator and denominator. Then when you subs in actual times, check whether the denominator is 0 before you divide by it.

[Mohammad Zainullah Khan]

Walter, you are a genius!

Did you manage to spot where I went wrong in the code above?

[Walter Roberson]

subs(solutions.x,{t1,t2},{T1 T2} )

[Mohammad Zainullah Khan]

That's prefect! Thanks.