Translating a Linear Regression from Matlab to Python
<p>I tried to translate a piece of code from Matlab to Python and I'm running into some errors:</p>

<p>Matlab:</p>

<pre><code>function [beta] = linear_regression_train(traindata)
y = traindata(:,1); %output
ind2 = find(y == 2);
ind3 = find(y == 3);
y(ind2) = -1;
y(ind3) = 1;
X = traindata(:,2:257); %X matrix,with size of 1389x256
beta = inv(X'*X)*X'*y;
</code></pre>

<p>Python:</p>

<pre><code>def linear_regression_train(traindata):
        y = traindata[:,0] # This is the output
        ind2 = (labels==2).nonzero()
        ind3 = (labels==3).nonzero()
        y[ind2] = -1
        y[ind3] = 1
        X = traindata[ : , 1:256]
        X_T = numpy.transpose(X)
        beta = inv(X_T*X)*X_T*y
        return beta
</code></pre>

<p>I am receiving an error: operands could not be broadcast together with shapes (257,0,1389) (1389,0,257) on the line where beta is calculated.</p>

<p>Any help is appreciated!</p>

<p>Thanks!</p>

<p>The problem is that you are working with numpy arrays, not matrices as in MATLAB.  Matrices, by default, do matrix mathematical operations.  So <code>X*Y</code> does a matrix multiplication of <code>X</code> and <code>Y</code>.  With arrays, however, the default is to use element-by-element operations.  So <code>X*Y</code> multiplies each corresponding element of <code>X</code> and <code>Y</code>.  This is the equivalent of MATLAB's <code>.*</code> operation.</p>

<p>But just like how MATLAB's matrices can do element-by-element operations, Numpy's arrays can do matrix multiplication.  So what you need to do is use numpy's matrix multiplication instead of its element-by-element multiplication.  For Python 3.5 or higher (which is the version you should be using for this sort of work), that is just the <code>@</code> operator.  So your line becomes:</p>

<pre><code>beta = inv(X_T @ X) @ X_T @ y
</code></pre>

<p>Or, better yet, you can use the simpler <code>.T</code> transpose, which is the same as <code>np.transpose</code> but much more concise (you can get rid of the `np.transpose line entirely):</p>

<pre><code>beta = inv(X.T @ X) @ X.T @ y
</code></pre>

<p>For Python 3.4 or earlier, you will need to use <code>np.dot</code> since those versions of python don't have the <code>@</code> matrix multiplication operator:</p>

<pre><code>beta = np.dot(np.dot(inv(np.dot(X.T, X)), X.T), y)
</code></pre>

<p>Numpy has a matrix object that uses matrix operations by default like the MATLAB matrix.  <em>Do not use it!</em>  It is slow, poorly-supported, and almost never what you really want.  The Python community has standardized around arrays, so use those.</p>

<p>There may also be some issues with the dimensions of <code>traindata</code>.  For this to work properly then <code>traindata.ndim</code> should be equal to <code>3</code>.  In order for <code>y</code> and <code>X</code> to be 2D, <code>traindata</code> should be <code>3D</code>.  </p>

<p>This could be an issue if <code>traindata</code> is 2D and you want <code>y</code> to be MATLAB-style "vector" (what MATLAB calls "vectors" aren't really vectors).  In numpy, using a single index like <code>traindata[:, 0]</code> reduces the number of dimensions, while taking a slice like <code>traindata[:, :1]</code> doesn't.  So to keep <code>y</code> 2D when <code>traindata</code> is 2D, just do a length-1 slice, <code>traindata[:, :1]</code>.  This is exactly the same values, but this keeps the same number of dimensions as <code>traindata</code>.</p>

<p><strong>Notes</strong>:  Your code can be significantly simplified using logical indexing:</p>

<pre><code>def linear_regression_train(traindata):
    y = traindata[:, 0] # This is the output
    y[labels == 2] = -1
    y[labels == 3] = 1
    X = traindata[:, 1:257]
    return inv(X.T @ X) @ X.T @ y
    return beta
</code></pre>

<p>Also, your slice is wrong when defining <code>X</code>.  Python slicing excludes the last value, so to get a 256 long slice you need to do <code>1:257</code>, as I did above.</p>

<p>Finally, please keep in mind that modifications to arrays inside functions carry over outside the functions, and indexing does not make a copy.  So your changes to <code>y</code> (setting some values to <code>1</code> and others to <code>-1</code>), will affect <code>traindata</code> outside of your function.  If you want to avoid that, you need to make a copy before you make your changes:</p>

<pre><code>y = traindata[:, 0].copy()
</code></pre>

python matlab