<div dir="ltr">I would certainly use einsum. It is almost perfect for these use cases, e.g.,<div>np.einsum('ki,kij,kj->k', A, inv(B), A)</div></div><br><div class="gmail_quote"><div dir="ltr">On Thu, Oct 26, 2017 at 12:38 PM Charles R Harris <<a href="mailto:charlesr.harris@gmail.com">charlesr.harris@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Oct 26, 2017 at 12:11 PM, Daniele Nicolodi <span dir="ltr"><<a href="mailto:daniele@grinta.net" target="_blank">daniele@grinta.net</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello,<br>
<br>
is there a better way to write the dot product between a stack of<br>
matrices?  In my case I need to compute<br>
<br>
y = A.T @ inv(B) @ A<br>
<br>
with A a 3x1 matrix and B a 3x3 matrix, N times, with N in the few<br>
hundred thousands range.  I thus "vectorize" the thing using stack of<br>
matrices, so that A is a Nx3x1 matrix and B is Nx3x3 and I can write:<br>
<br>
y = np.matmul(np.transpose(A, (0, 2, 1)), np.matmul(inv(B), A))<br>
<br>
which I guess could be also written (in Python 3.6 and later):<br>
<br>
y = np.transpose(A, (0, 2, 1)) @ inv(B) @ A<br>
<br>
and I obtain a Nx1x1 y matrix which I can collapse to the vector I need<br>
with np.squeeze().<br>
<br>
However, the need for the second argument of np.transpose() seems odd to<br>
me, because all other functions handle transparently the matrix stacking.<br>
<br>
Am I missing something?  Is there a more natural matrix arrangement that<br>
I could use to obtain the same results more naturally?</blockquote><div><br></div></div></div></div><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>There has been discussion of adding a operator for transposing the matrices in a stack, but no resolution at this point. However, if you have a stack of vectors (not matrices) you can turn then into transposed matrices like `A[..., None, :]`, so `A[..., None, :] @ inv(B) @ A[..., None]`  and then squeeze.</div><div><br></div><div>Another option is to use einsum.</div><div><br></div><div>Chuck</div></div></div></div>
_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@python.org" target="_blank">NumPy-Discussion@python.org</a><br>
<a href="https://mail.python.org/mailman/listinfo/numpy-discussion" rel="noreferrer" target="_blank">https://mail.python.org/mailman/listinfo/numpy-discussion</a><br>
</blockquote></div>