Dot product result has wrong dimension

Question

I'm doing a small project and I've encountered a problem while using Mathematica. I should mention that I'm kinda new to this software, so maybe it's just something I haven't grasped yet.

Anyway, I have a 9x9 matrix P and a 9D vector q0Dot.

Evaluating

projected = P.q0Dot

ended up taking a very long computational time. My laptop is quite powerful, but after I waited for like 10-15 minutes, I aborted the evaluation. I was kinda expecting slow evaluation, since the computational time of products is not linear.

So I reverted to calculate one element at a time, like:

a = P[[1, All]] . q0Dot

and so on (just one is enough for the sake of this post).

Now, the problem I've encountered is that a (and other elements of the result vector) is not a 9D vector, as it should be. is not a scalar.

I don't think it's related to the fact that Mathematica doesn't treat lists as row or column vectors. Correct me if I'm wrong.

To show what happened:

Dimensions[P]
Dimensions[q0Dot]
Dimensions[a]

{9, 9}
{9}
{8}

Note that almost every element of P cointains sine and cosine expressions and is fairly complicated.

So I was wondering what the problem is and made a couple of theories, but have found out nothing so far.

• maybe one element (maybe the last one) was null, so it got ignored at the end, but that's not the case;
• maybe the computational time is exceeding the maximum allowed and it stops computing without a warning, but I don't really know how to check that;
• maybe there's some Mathematica behaviour I don't know about

I'd be really really grateful if any of you could point me out at the solution!

EDIT: I'm adding my code, which I didn't do at first because I thought it would have been too detailed for the post itself. Also, it's making use of an external library, which I'm adding too.
My code + library - Google Drive

Now: it's an Inverse Kinematics problem, dealing with an industrial robotic arm. The matrix P (along with the others) contains only sines and cosines, but quite an amount, since it has 9 variables.

Unless you're doing something considerably more complicated than just having a dot product of integers or Reals, It should be super fast to calculate

In fact, I've tried that and didn't encounter any problems, as morbo stated.

Now I'm a bit confused. P . q0Dot should be a 9D vector, since it's a 9x9 matrix multiplied by a 9D (column) vector. As all of you stated, multiplying one row of P, which is a 9D (row) vector, by q0Dot should result in a scalar.
_{Well, blame my stupidty for this one, since I made a gaffe (which I'm going to highlight): maybe I didn't focus enough, since this should be obvious, duh.}
So why isn't it? Now I'm concerned about two things: one, is why that doesn't result in a scalar; two, what kind of computation Mathematica did to achieve that 8D result. Which may be the same problem in my case.

I think this is everything you asked for, to improve my question. If not, I'll surely add some more!

Answer 1

It'd be great if you posted/would show your code and matrices to see exactly what you're doing. It's a good habit to post your attempt (and your data) for us to play with.

Unless you're doing something considerably more complicated than just having a dot product of integers or Reals, It should be super fast to calculate,

matA = Table[RandomReal[], {i, 1, 9}, {j, 1, 9}];
vectorB = Table[RandomReal[], {i, 1, 9}];
vectorC = MatrixForm[Dot[matA, vectorB]]

It took me 0.000014 to calculate...even a 20k x 20k matrix took less than a second.

Your second bit of code is giving you a scalar because you are doing a scalar product of the first column of your 9x9 matrix with a 9 element vector, thus giving you a scalar product.

With Table you can populate the scalar products into the solution:

test = Table[matA[[i, All]] . matB, {i, 1, 9}] // MatrixForm

This gives the same result, which can be verified so:

vectorC == test

(*True*)

So, I'm guessing you're code is a bit wonky. Please give your working (or none working code) and I will update accordingly.

Addendum

I just saw that you're also dealing with some trig functions, without knowing which they are, I made a completely symbolic example

params = {x -> RandomReal[], y -> RandomReal[], z -> RandomReal[]};
AbsoluteTiming[
 var = {Cos[x], Cos[x] + Sin[y], Tanh[z], ArcSin[x], E^z, 
   ArcSinh[x] Cos[y]};
 smatA = Table[var[[RandomInteger[{1, 6}]]], {i, 1, 9}, {j, 1, 9}];
 svectorB = Table[var[[RandomInteger[{1, 6}]]], {i, 1, 9}];
 svectorC = MatrixForm[Dot[smatA, svectorB]];
 test = Table[smatA[[i, All]].svectorB, {i, 1, 9}] // MatrixForm;]
svectorC == test /. params

(*{0.000805, Null}*)
(*True*)

This also took via AbsoluteTiming[] 0.000805 seconds to generate randomly all equations, params, and solve.

; If you didn't know is shorthand for no output.

Answer 2

Edit

Since morbo has corrected his code, I no longer have any criticism of it to make, and simply offer this as alternative way of coding which is more concise and a little faster.

AbsoluteTiming[
  exprs = {Cos[x], Cos[x] + Sin[y], Tanh[z], ArcSin[x], E^z, ArcSinh[x] Cos[y]}; 
  matA = RandomChoice[exprs, {9, 9}];
  vecB = RandomChoice[exprs, 9];
  dot1 = matA.vecB;
  dot2 = Table[matA[[i]].vecB, {i, 9}];
  dot1 == dot2]

{0.000688, True}

Answer 3

Note:

Head[P[[1]].q0Dot]
(*  Plus  *)

And:

ecnt = 0;
Replace[P[[1]].q0Dot, e_ :> ToString[Head[e]][++ecnt], 1]
(*
"Times"[1] + "Times"[2] + "Times"[3] + "Times"[4] + "Times"[5] + 
 "Times"[6] + "Times"[7] + "Times"[8]
*)

What's going on is that Dimensions[] is returning the dimensions of the expression with the head Plus, which expression consists of 8 terms being added. Dimensions works with any head, not just List.

ArrayQ[P]
ArrayQ[q0Dot]
ArrayQ[a]
(*
  True
  True
  False
*)

Finally consider:

ByteCount[P] // AbsoluteTiming
(*  {10.5235, 19394763312}  *)

That's 19GB, which is more than my system can handle effectively. It also takes 10.5 sec for a simple operation. Currently, many of the subexpressions appear to be shared (since the MemoryInUse is a couple hundred MB). Almost any operation that causes P to be traversed (such as Dot) would take a long time, and if subexpressions need copying, it would require swapping. The kernel might even run out of memory and crash.

Maybe if you have 64+GB of RAM, things will work better. Otherwise, see if you can deal more effectively with P. One culprit is taking the Inverse[] of a symbolic matrix. It's unlikely that one would be able to analyze such large symbolic expressions, so if you're ultimately going to do things numerically, then wait to construct P when all parameters are numeric. (Then you'll probably want to use LinearSolve[] instead of Inverse[].) It should be much, much faster.

Answer 4

Here's a fix

Although I agree with @MichaelE2 that inverting a symbolic matrix is asking for trouble!

First, in your W matrix, write 1/2 instead of 0.5 and some odd numerical behaviors downstream are averted. When it is time to simplify things, the algorithms can do a little better.

Second, it doesn't hurt to keep simplifying your code as you go, so for example using FullSimplify early and simplifying JRNonInv helps keep the bloat down when you invert the matrix.

JTransposed = Transpose[J];
JSimplified = J // FullSimplify;  (* was just a simplify here *)
JTSimplified = Transpose@JSimplified ;   (* no need to simplify JTransposed *)  

JRNonInv = (JSimplified.WInverse.JTSimplified) // Simplify; (* added a simplify *)

To your point on the code not working, if you replace q0Dot with a "proper" vector...

qvec = Transpose@{q0Dot};
qvec // Dimensions

{9,1}

Your code now works, at least in creating the correct sized matrices.

Dimensions[P]
Dimensions[qvec]
Dimensions[a]

{9, 9}

{9, 1}

{1}

You can also skip the intermediate variables and go

 projected = P.qvec
 Dimensions@projected

{9,1}

I firmly believe how Mathematica allows you to define a vector as a simple list when dealing with linear algebra is the Devil's work, and has caused many a heartbreak for me.

EDIT

@ slim71, here's an example of what I mean. Define a simple "vector"

vsimp = {1,1,1};
vsimp//TeXForm

$$ \{1,1,1\} $$

If you just do MatrixForm it looks like a vector, but it's not really one. Now make a full-blown column vector, which is a $3 \times 1$ matrix.

vector = Transpose@{{1, 1, 1}};
vector//TeXForm

$$ \left( \begin{array}{c} 1 \\ 1 \\ 1 \\ \end{array} \right) $$

Create a simple matrix to play with

mat = Table[i,{i,3},{3}];
mat//TeXForm

$$ \left( \begin{array}{ccc} 1 & 1 & 1 \\ 2 & 2 & 2 \\ 3 & 3 & 3 \\ \end{array} \right) $$

Now in theory we can do mat.vector but not vector.mat since the dimensions would be incompatible in the second case.

mat.vector

$$ \left( \begin{array}{c} 3 \\ 6 \\ 9 \\ \end{array} \right) $$

The command vector.mat fails, complaining about incompatible shapes. Give it a try.

In the other case, we can execute both vsimp.mat and mat.vsimp, depending on Mathematica to interpret vsimp either as a row vector or column vector given the context. You can easily get burned with a mistake here. Writing the vector as a matrix helps to prevent these mistakes, and it appears to fix the error in your code, which ought to work correctly. Not sure why it doesn't, other than that you've got a monstrous term for the inverse of the matrix.

Finally, you can't just do Transpose@vsimp since you transpose matrices, and vsimp is just a list.

Transpose@vsimple
(* Transpose::nmtx: The first two levels of {1,1,1} cannot be transposed. *)

Instead you are constructing a matrix to represent your column vector.