Suppose you have n square matrices A1,...,An. Is there anyway to multiply these matrices in a neat way? As far as I know dot in numpy accepts only two arguments. One obvious way is to define a function to call itself and get the result. Is there any better way to get it done?
6 Answers
This might be a relatively recent feature, but I like:
A.dot(B).dot(C)
or if you had a long chain you could do:
reduce(numpy.dot, [A1, A2, ..., An])
Update:
There is more info about reduce here. Here is an example that might help.
>>> A = [np.random.random((5, 5)) for i in xrange(4)]
>>> product1 = A[0].dot(A[1]).dot(A[2]).dot(A[3])
>>> product2 = reduce(numpy.dot, A)
>>> numpy.all(product1 == product2)
True
Update 2016:
As of python 3.5, there is a new matrix_multiply symbol, @
:
R = A @ B @ C
-
Thanks for the response. The first option works fine; but the second one doesnt; or at least I couldn't make it work. Can you please elaborate it a bit more or maybe give an example? Thanks a lot– NNsrCommented Aug 7, 2012 at 5:41
-
7I run into this all the time and ended up writing a helper function. Wish this was part of NumPy:
def xdot(*args): return reduce(np.dot, args)
– rd11Commented Jul 3, 2014 at 9:57 -
1
-
4In python3 reduce was moved to functools,
from functools import reduce
.– Bi RicoCommented May 15, 2020 at 19:28 -
1Can you edit to put
R = A @ B @ C
at the top of the answer (I cannot edit because suggested queue is full. Then duckduckgo will preview it when this question is searched. Right now, it shows
A.dot(B).dot(C)` Commented Feb 2, 2022 at 12:50
Resurrecting an old question with an update:
As of November 13, 2014 there is now a np.linalg.multi_dot
function which does exactly what you want. It also has the benefit of optimizing call order, though that isn't necessary in your case.
Note that this available starting with numpy version 1.10.
-
2This should be the go-to answer for anyone reaching here in the future. Commented Apr 24, 2020 at 20:06
-
1Is there a speed difference between
@
andnp.linalg.multi_dot
?– kanso37Commented Aug 23, 2020 at 15:43 -
3@kanso37 I created a list of arrays as shown above using
A_list = [np.random.random(100, 100) for i in range(3)]
run a simple test using%timeit np.linalg.multi_dot(A_list)
vs%timeit A_list[0] @ A_list[1] @ A_list[2]
. It seems that the second method outperforms the first one (100 us vs 85 us on my machine), but of course I cannot tell whether this is true in general. I also wonder how the second method would be generalized using the list recursively as the first one. Commented Nov 17, 2020 at 17:25
If you compute all the matrices a priori then you should use an optimization scheme for matrix chain multiplication. See this Wikipedia article.
-
2Thanks for your comment; but I dont think for square matrices it matters. Right?– NNsrCommented Aug 7, 2012 at 2:35
-
@Nikandish: Correct. I missed that part in your original answer. Commented Aug 8, 2012 at 14:57
Another way to achieve this would be using einsum
, which implements the Einstein summation convention for NumPy.
To very briefly explain this convention with respect to this problem: When you write down your multiple matrix product as one big sum of products, you get something like:
P_im = sum_j sum_k sum_l A1_ij A2_jk A3_kl A4_lm
where P
is the result of your product and A1
, A2
, A3
, and A4
are the input matrices. Note that you sum over exactly those indices that appear twice in the summand, namely j
, k
, and l
. As a sum with this property often appears in physics, vector calculus, and probably some other fields, there is a NumPy tool for it, namely einsum
.
In the above example, you can use it to calculate your matrix product as follows:
P = np.einsum( "ij,jk,kl,lm", A1, A2, A3, A4 )
Here, the first argument tells the function which indices to apply to the argument matrices and then all doubly appearing indices are summed over, yielding the desired result.
Note that the computational efficiency depends on several factors (so you are probably best off with just testing it):
A_list = [np.random.randn(100, 100) for i in xrange(10)]
B = np.eye(A_list[0].shape[0])
for A in A_list:
B = np.dot(B, A)
C = reduce(np.dot, A_list)
assert(B == C)
This works in VS Code for two matrices
import numpy as np #
def matrix_multiply(matrix1, matrix2):
print(f"Matrix A:\n {A}\n")#Print the Matrix contents
print(f"Matrix B:\n {B}\n")
if A.shape[1] == B.shape[0]:#Check if matrices can be multiplied
C = np.matmul(A,B) #Use matmul to multiply the matrices
return C #Return the resulting matrix
else:
return "Sorry, cannot multiply A and B."#Error catching
# Use np to generate dataset
np.random.seed(27)
A = np.random.randint(1,10,size = (5,4))
B = np.random.randint(1,10,size = (4,2))
# Testing the function
result= matrix_multiply(A,B) #Call matrix_multiply to find answer
print(result)
#References:
#https://geekflare.com/multiply-matrices-in-python/#geekflare-toc-use-python-nested-list-comprehension-to-multiply-matrices
#https://www.anaconda.com/download
#Launch VS Code from Anaconda