How do we know the relation between times and the number of arithmetic floating point ?

I am trying to find out if solution of linear systems of equations with a triangular matrix by the backslash operation in Matlab requires about n^2 or n^3 arithmetic floating point operations, but I do know how to find the follwing:

How to compute the average time for the computation x=L\b using the tic-toc?

how to plot a graph of the average times versus n using semilogy?

what the exponent of n is?

Can you please explain it to me? Thank you in advance.

n=100;
L=tril(randn(n));
new=ones(1,n);
L(1,:)=new
b=ones(n,1);
x=L\b;

[Rik]

The tic,toc construct should not be used for accurate performance estimation, use timeit instead. If it warns the operation is too fast, you can put the operation in a loop.

To plot the trend you will need to store the time in an array.

[Omar B.]

I am new to Matlab, and do not know how to do that.

[Rik]

If you have trouble with Matlab basics you may consider doing the Onramp tutorial (which is provided for free by Mathworks).

[Omar B.]

Thank you for your advice.

Answer 1

[Walter Roberson]

Generally speaking, measure the time for a calculation on a matrix m x m, then for one a*m x a*m, then b*m x b*m. If the calculation is linear in m, then the time for the first should be m^2*C1+C2 for constants C1 and C2, and the time for the middle should be a^2*m^2*C1+C2 and that is enough information to solve for C1 and C2. Use those to make a prediction for b*m. Was the third time close to that? If so then the time is approximately linear in m^2. Or quadratic in m, depending on how you want to look at it.

Now if the times are not predictive for b*m then change your model ffrom ^2 to ^n for unknown n, and write out the math for m, a*m, and b*m. Three unknowns, three equations, and you can proceed to estimate n. You can cross-check with a fourth test.

[Omar B.]

I did not get waht you mean. I am new to Matlab. How can we show this via poltting the time?

[Walter Roberson]

m = 50;
a = 5;
b = 20;
tic;
m1 = sprand(1*m, 1*m, 0.01);
t1 = toc;
tic;
m2 = sprand(a*m, a*m, 0.01);
t2 = toc;
tic;
m3 = sprand(b*m, b*m, 0.01);
t3 = toc;
plot([1, a, b], [t1, t2, t3])

semilogy([1, a, b], [t1, t2, t3])

p = polyfit([1, a, b], [t1, t2, t3], 2)

p = 1×3

0.0003 -0.0068 0.0465

[Omar B.]

Thank you, but I do not get what the relation between what you did and solving system of linear equation.

[Walter Roberson]

polyfit(x, y, 2) effectively solves a series of linear equations that are of the form

p2*x.^2 + p1*x + p0 = y

These are linear equations because x and y are known. In this case,

p2*1^2 + p1*1 + p0 = t1
p2*5^2 + p1*5 + p0 = t2
p2*20^2 + p1*20 + p0 = t3

which is

  1*p2 +  1*p1 + 1*p0 = t1
 25*p2 +  5*p1 + 1*p0 = t2
400*p2 + 20*p1 + 1*p0 = t3

which could also be written as

[  1  1 1  * [p2; = [t1;
  25  5 1     p1;    t2;
 400 20 1]    p0]    t3]

or

A = [1 1 1; 25 5 1; 400 20 1];
x = [p2; p1; p0]
b = [t1; t2; t3]

and A*x = b, which leads you to

A\b

which is

[1 1 1; 25 5 1; 400 20 1] \ [t1; t2; t3]

You could certainly argue with my choice of using [1, a, b] as the x coordinates for the fit. You could, for example, instead use [m, a*m, b*m], or [m^2, a^2*m^2, b^2*m^2] .

Mathematically, the [1,a,b] and [m,a*m,b*m] cases are equivalent to a constant of proportionality. Both of them treat the situation as if the "size" of the operation depended on the number of number of rows (or columns) of the square matrix, L in your code -- as if perhaps the cost of doing the \ operator depended on your n in randn(n).

Now think about

p3*x.^3 + p2*x.^2 + p1*x + p0 = y

and this time you would probably want to use [m, m*a, m*b, m*c] and [t1, t2, t3, t4] as coordinates for size tests of m x m, (a*m) x (a*m), (b*m) x (b*m), (c*m) x (c*m) . Create the simultaneous equations (all of the x and y values are known so you are arriving at constants times p3, p2, p1, p0), solve... and you get a trial polynomial that tries to directly relate the cost of the operation you are testing to the "order" of square you are testing with.

If you are wondering what sprand() has to do with this: sprand() is just some operation that I was showing you could plot the timing for, since what you asked for an illustration of was the plotting part.

[Omar B.]

Thank you.