lsmr not being used on multidimensional regression

[grudloff]

In PR #9 for multidimensional output. scipy's lsmr only accepts b as a vector, so in this scenario, it raises a ValueError and it is completely bypassed and the more expensive pinv is used.

I think this should be addressed. I am thinking of tackling this by using numpy's lstsq that supports 2-dimensional 'b', or the other alternative would be to use lsmr on each slice of 'b'. Something like this:

try:
    if len(shape)>1:
        z = np.linalg.lstsq(D.T, t, rcond=None)[0]
    else:
        z = linalg.lsmr(D.T, t)[0]
except:
    z = np.linalg.pinv(D).T @ t

or if we stick to lsmr, something like this:

try:
    if len(shape)>1:
        z = np.column_stack([linalg.lsmr(D.T, t_slice)[0] for t_slice in t.T])
    else:
        z = linalg.lsmr(D.T, t)[0]
except:
    z = np.linalg.pinv(D).T @ t