differences between 2 identical instructions

Hi,

I've got a problem with Simulink. I've created a Simulink model, which takes a matrix (400 x 10). The model transpose that matrix, after does a matrix moltiplication with a vector (400,1). There is a difference between this output and the output of the matlab file, which does the following operations:

vect_out= U' * vect_in;

the difference is 10^-12, but it causes problems in the following steps.

There is a difference in the results also if in matlab I write these instructions:

A= U'; vect_out1= A* vect_in;

vect_out-vect_out1 is not a zeros vector, the difference is 10^-12. U'*vect_in give me the correct results at the end of the matlab program (all the matlab program, not only the instructions that I've written here) . When I break that operation in two operations (Simulink does that) the result (A= U'; vect_out1= A* vect_in) is different.

how can I implement U'*vect in simulink in only one operation? What are the causes of that difference in the results? How can I resolve that problem?

Answer 1

[Jan]

How do you know which of the results is "correct"? The matrix multiplication is prone to rounding effects. A trivial example is:

a = [1, 1, -1]
b = [1e17, 1, 1e17]
a * b'

Now imagine than the order of calculations is different for:

x1 = A' * y;

and

B  = A';
x2 = B * y;

because in the first case a BLAS function calculates this without an explicit transposition. But neither x1 nor x2 are more correct, but both are affected by rounding errors.

If an algorithm critically depends on such rounding effects, it is called "numerically instable" and a scientifically hard conclusion is, that the results are random.

[Teja Muppirala]

To expand on what Jan is saying, in the numerical world, there are all sorts of weird things that ideally should not happen, but actually they do due because order of operations changes results.

Example 1 (this pertains to your problem):

A = rand(100);
x = rand(100,1);
isequal( A'*x,  (A'+0)*x)

ans =

       0

Example 2.

r1 = randn(1e6,1);
r2 = randn(1e6,1);

[r1'*r2] - [flipud(r1)' * flipud(r2)]

ans =

     5.9117e-12

Example 3.

A = randn(1000);
isequal(sum(A(:)), sum(sum(A)))

ans =

       0

This list goes on and on. To answer your specific question, I don't know how you would implement a "transpose multiply" in one Simulink operation that will give you the same answer as U'*v in MATLAB. But if your algorithm breaks down because of simple numerical errors like this, you may need come up with a more robust way of solving your problem.

[Jan]

@Teja: Exactly. In addition, sometimes the algorithm is numerically instable and an improvement solves the problem. But sometimes the problem is instable itself, e.g. a weather forecast for the next 20 years, a falling snow flake, a feeback controlled fast moving multi-body system like a running robot or cat. In general the state of such "instable" systems is influenced by external events (cosmic rays, falling ricebags in China, butterflies in Brasilia, depressive brokers in the Wallstreet and messy carry bits in the floating point units of processors). Therefore a scientifically meaningful simulation requires an analysis of the sensitivity of the results to small variations of the initial values and parameters. I frequently have to grin, when students simulate the flight of a spaceship back from Mars to Earth with the re-entrance in the atmosphere and are convinced that they have determined the ground position with an accuracy of 1e-8 meter. The principal limitations get immediately clear, when I ask for a validation of the simulation by measuring the final position in the real object. Then no student assumes, that a spaceship can be located with such an accuracy, because this is even physically meaningless due to the spatial dimensions of the object.

So a simulation has to consider 3 effects: the mathematical model, the implementation of the algorithm and the properties of the physical/real world object. Do not see the limitations of the linear algebra as an isolated source of problems.

Or in other words: "All sorts of weird things" happen in the numerical world and in the real world also. Therefore the definition of "weird" is very relative and should be replaced by "normal".