[Numpy-discussion] matrix default to column vector?

[Charles R Harris]

On Sun, Jun 14, 2009 at 10:20 AM, Tom K.<tpk at kraussfamily.org> wrote:
>
>
>
> jseabold wrote:
>>
>> On Mon, Jun 8, 2009 at 3:33 PM, Robert Kern<robert.kern at gmail.com> wrote:
>>> On Mon, Jun 8, 2009 at 14:10, Alan G Isaac<aisaac at american.edu> wrote:
>>>>>> Going back to Alan Isaac's example:
>>>>>> 1)  beta = (X.T*X).I * X.T * Y
>>>>>> 2)  beta = np.dot(np.dot(la.inv(np.dot(X.T,X)),X.T),Y)
>>>>
>>>>
>>>> Robert Kern wrote:
>>>>> 4) beta = la.lstsq(X, Y)[0]
>>>>>
>>>>> I really hate that example.
>>
>
> I propose the following alternative for discussion:
>
>  U, s, Vh = np.linalg.svd(X, full_matrices=False)
>  1)  beta = Vh.T*np.asmatrix(np.diag(1/s))*U.T*Y
>  2)  beta = np.dot(Vh.T, np.dot(np.diag(1/s), np.dot(U.T, Y)))
>
> 1) is with X and Y starting out as matrices, 2) is with .
>.  Sigh.
>

The problem is that I left the diagonal returned by svd as an array
rather than a matrix for backward compatibility. Diag returns the
diagonal when presented with a 2d matrix (array).  I went back and
forth on that, but not doing so would have required everyone to change
their code to use diagflat instead of diag. I do note that diag
doesn't preserve matrices and that looks like a bug to me.

This is also an argument against over-overloaded functions such as
diag. Such functions are one of the blemishes of MatLab, IMHO. OTOH,
there is something of a shortage of short, meaningful names

Chuck