Sorry, the linearly independent thing is a bit of a read herring. The
following code is equivalent to using np.linalg.solve if X is (b, b) and
np.linalg.lstsq if X is (b, m) for the 1d case, but it also works when data
is (x, y, z, b).

invX = np.linalg.pinv(X)
result1 = np.dot(data, invX.T)
# Or
result2 = np.tensordot(invX, data, axes=(1, -1))
# if data is (x, y, z, b) result 2 will be (b, x, y, z) to get it to be (x,
y, z, b) we can use rollaxis
# notice that result2.T will make it (b, z, y, x)
result2 = np.rollaxis(result2, 0, result2.ndim)


On Thu, Aug 22, 2013 at 9:20 AM, Ariel Rokem <notifications@github.com>wrote:

> In dipy/reconst/dki.py:
>
> > +   is given by:
> > +
> > +    ..math ::
> > +
> > +        \hat{\beta} = (A' x A)^{-1} A' y
> > +
> > +    See also
> > +    --------
> > +    http://en.wikipedia.org/wiki/Ordinary_least_squares#Estimation
> > +    """
> > +
> > +    A = np.asarray(A)
> > +
> > +    X = np.matrix(A.copy())
> > +
> > +    return np.dot(linalg.pinv(np.dot(X.T, X)), X.T)
>
> You are correct. I wanted a way to perform the operation of linalg.solve
> for a matrix with dimensionality (x,y,z, b). If we use linalg.solve we
> would have to loop over voxels, which is what I am trying to avoid here.
> Does that even make sense though?
>
> Here - the rows will sometime be linearly dependent - in cases where
> measurements in two b values are done along the same direction. I haven't
> tested that yet, but that might be a problem (?).
>
> —
> Reply to this email directly or view it on GitHub<https://github.com/nipy/dipy/pull/233/files#r5928786>
> .
>