Paolo de Leva

Shatrughan, I do not know outerm and PLS, and I cannot completely understand your request, but I can try and explain what you can do with my functions. Allow me to start with the simplest case:

A = rand(7,1);
B = rand(9,1);
X = outer(A, B);

Here, X is a matrix, and its size is 7x9. It would make sense to write a function to generalize this to 3-D:

A = rand(7,1);
B = rand(9,1);
C = rand(10,1);
X = simple_outer3D(A, B, C);

where X is a 3-D array and its size is 7x9x10, which is the size of your expected result. This can be done using this code (where BAXFUN is a function that I posted on MATLAB Central File Exchange; File ID: #23084):

function X = simple_outer3D(A, B, C)
% A, B, and C must be column vectors
X = baxfun(@times, A, B, 0, 1); % Size: 7x9
X = baxfun(@times, X, C, 0, 2); % Size: 7x9x10

Warning: this code is completely valid as a generalization of an outer product only if A, B, and C do not contain complex numbers. See the help of function OUTER for details.

Notice that I got the result you expected, but the input matrices have not the same size as in your example. Conversely, input matrices with the same size as in your example would give a different result. Simple case:

A = rand(7,2);
B = rand(9,2);
X = outer(A, B);

Here, A and B are not two vectors, but two “arrays of vectors” (each containing two vectors). Hence, X is not a single outer product, but a 7x9x2 array containing two outer products. The first outer product is X(:,:,1) and the second is X(:,:,2). If we generalize to 3-D, we have:

A = rand(7,2);
B = rand(9,2);
C = rand(10,2);
X = outer3D(A, B, C);

In this case, X is supposed to contain two 3-D arrays. Hence, its size should be 7x9x10x2 (rather than 7x9x10). This can be done by defining outer3D as follows:

function X = outer3D(A, B, C)
A = reshape(A, [size(A,1),1,1,size(A,2)]); % Size: 7x1x1x2
B = reshape(B, [1,size(B,1),1,size(B,2)]); % Size: 1x9x1x2
C = reshape(C, [1,1,size(C)]); % Size: 1x1x10x2
X = baxfun(@times, A, B); % Size: 7x9x1x2
X = baxfun(@times, X, C); % Size: 7x9x10x2

Warning: this code is completely valid as a generalization of an outer product only if A, B, and C do not contain complex numbers. See the help of function OUTER for details.
Note: In this code, BAXFUN (binary array expansion) is equivalent to BSXFUN (binary singleton expansion).

This is even better:

"DCM - 3x3xN multidimensional direction cosine matrix which performs coordinate system rotations by pre-multiplying column vectors (V_rot = DCM * V, where V_rot is V resolved in the rotated coordinate system)."

Paolo de Leva

OK, John,

I GUESS THIS IS OUR COMMON BACKGROUND:

We agree that the rotation from a reference frame A (original) to another reference frame B (rotated) produces an apparent rotation in the opposite direction of a vector v which is fixed in A. This leads to two conventions to describe rotation. Let’s call them:
1) "FRAME ROTATION"
2) “VECTOR ROTATION" (or “vector space rotation”)
The most important point is that, when we describe a rotation from A to B, we must assume that v is fixed in A (see below). Also, we agree that SpinCalc uses the “FRAME ROTATION” convention, and of course we agree on Rodriguez's rotation formula and all other conversion formulas which are correctly implemented in SpinCalc.

Moreover, as discussed in the previous message, we agree that there are two conventions used to represent rotation with a DCM (pre- or post-multiplication). Your recently updated help text makes clear you are using pre-multiplication of a column vector by the DCM. IF we call R the DCM, and v the original column vector, then R is a matrix such that:

v_rot = R*v.

NOW, HERE'S WHAT WE DISAGREE ABOUT:
 
We do not agree about your help text. In your second last message, you explained that you implemented the “FRAME ROTATION” convention, but this is not what you wrote in your help text, which defines the DCM as:

"a multidimensional matrix which pre-multiplies a coordinate frame (column) vector to rotate it to the desired new frame."

Let’s call A the initial frame, and B “the desired new frame“.
It unequivocally follows from your text that, in the equation v_rot = R*v, you describe v as one of the versors of A (i.e., a standard basis vector; i.e., either [1,0,0], [0,1,0], or [0,0,1], in 3D). Actually, in this singular case R*v is the same as selecting one of the columns of R, and this description of v is quite questionable, in my opinion, as a (pre- or post-) multiplication by R is typically used to rotate any vector, not only the versors of A. But this is not the most important fault in your definition.

Since your help text describes the rotation of a "column vector", then you unequivocally defined R as a matrix which produces (i.e. represents) a “VECTOR ROTATION”, not a “FRAME ROTATION”. Moreover, the “FRAME ROTATION” convention assumes that v “is fixed in A” (see above). On the contrary, according to your definition, v unequivocally appears fixed in B (not in A).

Most importantly and unfortunately, R does not represent a “FRAME ROTATION” from A to B, but a “FRAME ROTATION” of equal magnitude in the opposite sense, from A to C. Namely, C does not coincide with B, as B rotates together with the vector. For instance, if B is rotated by 45° about z, then C is rotated by -45° about the same axis. Thus, we can rewrite the equation as follows:

vC = R*vA.
(vB = R*vA would be incorrect).

In short, your help text describes a “FRAME ROTATION” from A to B, which is the same as (not the opposite of) a the “VECTOR ROTATION” rotation produced by R (from v to v_rot, or from vA to vC, if you prefer). But unfortunately R neither produces nor represents the frame rotation from A to B.

However, your axis-angle representation describes the opposite rotation, the “FRAME ROTATION” from A to C, in which v is “fixed in A” and observed from C. This is the frame rotation produced and represented by R, but C is not even mentioned in your help text! In other words, your definition of rotation is inconsistent with your implementation, and this is higly misleading.

Paolo de Leva

Paolo

Shatrughan, what is the code you would write without having multiprod()?

As far as I can see, multiprod() accepts 2 input arrays while you want to do something with 3 arrays.

multiprod() does nothing which you cannot do with regular code, it just does it more elegantly and faster.

So, please write down an example code which you expect to get "beautified" or accelerated with the help of multiprod().

This is great -- it has been a tremendous improvement for my project.