59

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 6

88

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
6
  • 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
    – NNsr
    Commented Aug 7, 2012 at 5:41
  • 7
    I 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)
    – rd11
    Commented Jul 3, 2014 at 9:57
  • 1
    NameError: name 'reduce' is not defined
    – rsc05
    Commented Mar 2, 2020 at 13:57
  • 4
    In python3 reduce was moved to functools, from functools import reduce.
    – Bi Rico
    Commented May 15, 2020 at 19:28
  • 1
    Can 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
53

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.

3
  • 2
    This should be the go-to answer for anyone reaching here in the future.
    – SajidSalim
    Commented Apr 24, 2020 at 20:06
  • 1
    Is there a speed difference between @ and np.linalg.multi_dot?
    – kanso37
    Commented 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
5

If you compute all the matrices a priori then you should use an optimization scheme for matrix chain multiplication. See this Wikipedia article.

2
  • 2
    Thanks for your comment; but I dont think for square matrices it matters. Right?
    – NNsr
    Commented Aug 7, 2012 at 2:35
  • @Nikandish: Correct. I missed that part in your original answer. Commented Aug 8, 2012 at 14:57
4

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):

2
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)
0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.