Matrices addition in numpy
<p>I need to calculate the sum of a list of matrices, however, I can't use <code>np.sum</code>, even with <code>axis=0</code>, I don't know why. The current solution is a loop, but is there a better way for that? </p>

<pre><code>import numpy as np
SAMPLE_SIZES = [10, 100, 1000, 10000]
ITERATIONS = 1
MEAN = np.array([1, 1])
COVARIANCE = np.array([[1, 0.5], [0.5, 1]])
for sample_size in SAMPLE_SIZES:
    max = -1
    for i in range(ITERATIONS):
        xs = np.random.multivariate_normal(MEAN, COVARIANCE, size=sample_size)
        sigma = [[0, 0], [0, 0]]
        for x in xs:
            sigma += np.outer((x-MEAN), (x-MEAN)) / (sample_size-1)
</code></pre>

<p>In the code above, can I replace the last loop using some <code>numpy</code> function? I guess using a loop would be not efficient if the data is very large.</p>

<p>Read up about numpy <a href="https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html" rel="nofollow">broadcasting</a>.</p>

<pre><code>xs = np.random.multivariate_normal(MEAN, COVARIANCE, size=sample_size)
</code></pre>

<p><code>xs</code> now has shape <code>(sample_size, 2)</code>, which means you can just subtract <code>MEAN</code> directly. You now need to take the outer product between <code>xs - MEAN</code> and <code>xs - MEAN</code> while adding over the <code>sample_size</code> axis. This is best done using <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html" rel="nofollow"><code>np.einsum</code></a>:</p>

<pre><code>&gt;&gt;&gt; sigma = np.einsum('ij,ik-&gt;jk', xs - MEAN, xs - MEAN) / sample_size
&gt;&gt;&gt; sigma    
array([[ 1.00216043,  0.49549231],
       [ 0.49549231,  1.00004423]])
</code></pre>

<p>An alternative is to use broadcasting:</p>

<pre><code>&gt;&gt;&gt; sigma = np.sum((xs - MEAN)[:, :, np.newaxis]
                   * (xs - MEAN)[:, np.newaxis, :], axis=0) / sample_size
</code></pre>

<p>Though the broadcasting solution seems easier to understand, <code>np.einsum</code> will usually be <a href="http://stackoverflow.com/a/39482794/525169">more efficient</a> than broadcasting.</p>

<p><strong>Additional note:</strong> Note that I've divided by <code>sample_size</code>, and not by <code>sample_size - 1</code>. This is because for estimating the covariance matrix of a random variable with known mean, you need to divide by <code>sample_size</code>. Use <code>sample_size - 1</code> when you are estimating the mean from the same dataset as well, and using it in your covariance estimate. Otherwise your covariance estimate will be biased.</p>

python numpy