memory error while performing matrix multiplication
<p>As part of a project I'm working on, I need to calculate the mean squared error between <code>2m</code> vectors.</p>

<p>Basically I two matrices <code>x</code> and <code>xhat</code>, both are size <code>m</code> by <code>n</code> and the vectors I'm interested in are the rows of these vectors.</p>

<p>I calculate the MSE with this code</p>

<pre><code>def cost(x, xhat): #mean squared error between x the data and xhat the output of the machine
    return (1.0/(2 * m)) * np.trace(np.dot(x-xhat,(x-xhat).T))
</code></pre>

<p>It's working correctly, this formula is correct.</p>

<p>The problem is that in my specific case, my <code>m</code> and <code>n</code> are very large. specifically, <code>m = 60000</code> and <code>n = 785</code>. So when I run my code and it enters this function, I get a memory error.</p>

<p>Is there a better way to calculate the MSE? I'd rather avoid for loops and I lean heavily towards matrix multiplication, but matrix multiplication seems extremely wasteful here. Maybe something in numpy I'm not aware of?</p>

<p>The expression <code>np.dot(x-xhat,(x-xhat).T)</code> creates an array with shape (m, m).  You say m is 60000, so that array is almost 29 gigabytes.</p>

<p>You take the trace of the array, which is just the sum of the diagonal elements, so most of that huge array is unused.  If you look carefully at <code>np.trace(np.dot(x-xhat,(x-xhat).T))</code>, you'll see that it is just the sum of the squares of all the elements of <code>x - xhat</code>.  So a simpler way to compute <code>np.trace(np.dot(x-xhat,(x-xhat).T))</code> that doesn't require the huge intermediate array is <code>((x - xhat)**2).sum()</code>.  For example,</p>

<pre><code>In [44]: x
Out[44]: 
array([[ 0.87167186,  0.96838389,  0.72545457],
       [ 0.05803253,  0.57355625,  0.12732163],
       [ 0.00874702,  0.01555692,  0.76742386],
       [ 0.4130838 ,  0.89307633,  0.49532327],
       [ 0.15929044,  0.27025289,  0.75999848]])

In [45]: xhat
Out[45]: 
array([[ 0.20825392,  0.63991699,  0.28896932],
       [ 0.67658621,  0.64919721,  0.31624655],
       [ 0.39460861,  0.33057769,  0.24542263],
       [ 0.10694332,  0.28030777,  0.53177585],
       [ 0.21066692,  0.53096774,  0.65551612]])

In [46]: np.trace(np.dot(x-xhat,(x-xhat).T))
Out[46]: 2.2352330441581061

In [47]: ((x - xhat)**2).sum()
Out[47]: 2.2352330441581061
</code></pre>

<p>For more ideas about computing the MSE, see the <a href="http://stackoverflow.com/questions/16774849/mean-squared-error-in-numpy">link provided by user1984065</a> in a comment.</p>

python numpy math memory matrix