Estimation for many voxels at the same time

We often want to fit the same design to many different voxels.

Let’s make a design with a linear trend and a constant term:

import numpy as np
import matplotlib.pyplot as plt
# Print arrays to 2 decimal places
np.set_printoptions(precision=2)

X = np.ones((12, 2))
X[:, 0] = np.linspace(-1, 1, 12)
plt.imshow(X, cmap='gray')

<matplotlib.image.AxesImage at 0x7f53c1ccb1c0>

To fit this design to any data, we take the pseudo-inverse:

import numpy.linalg as npl
piX = npl.pinv(X)
piX

array([[-0.21, -0.17, -0.13, -0.1 , -0.06, -0.02,  0.02,  0.06,  0.1 ,
         0.13,  0.17,  0.21],
       [ 0.08,  0.08,  0.08,  0.08,  0.08,  0.08,  0.08,  0.08,  0.08,
         0.08,  0.08,  0.08]])

Notice the shape of the pseudo-inverse:

piX.shape

(2, 12)

Now let’s make some data to fit to. We will draw some samples from the standard normal distribution.

# Random number generator.
rng = np.random.default_rng()
# 12 random numbers from normal distribution, mean 0, std 1.
y_0 = rng.normal(size=12)
y_0

array([-0.91,  0.7 , -0.34, -1.56, -0.39,  0.03,  0.4 , -0.4 ,  0.19,
        1.32,  1.23, -0.24])

beta_0 = piX @ y_0
beta_0

array([0.63, 0.  ])

We can fit this same design to another set of data, using our already-calculated pseudo-inverse.

y_1 = rng.normal(size=12)
y_1

array([ 0.55, -0.1 , -0.52, -0.36, -0.74, -0.21,  0.79,  0.26, -2.12,
        1.15,  0.83, -0.98])

beta_1 = piX @ y_1
beta_1

array([-0.03, -0.12])

And another!:

y_2 = rng.normal(size=12)
beta_2 = piX @ y_2
beta_2

array([-0.47, -0.09])

Now the trick. Because of the way that matrix multiplication works, we can fit to these three sets of data with the one single matrix multiply:

# Stack the data vectors into columns in a 2D array.
Y = np.vstack((y_0, y_1, y_2)).T
Y

array([[-0.91,  0.55,  1.61],
       [ 0.7 , -0.1 ,  0.18],
       [-0.34, -0.52, -1.  ],
       [-1.56, -0.36, -0.35],
       [-0.39, -0.74, -0.66],
       [ 0.03, -0.21, -0.08],
       [ 0.4 ,  0.79,  1.47],
       [-0.4 ,  0.26, -0.71],
       [ 0.19, -2.12,  0.5 ],
       [ 1.32,  1.15, -0.95],
       [ 1.23,  0.83, -0.9 ],
       [-0.24, -0.98, -0.24]])

betas = piX @ Y
betas

array([[ 0.63, -0.03, -0.47],
       [ 0.  , -0.12, -0.09]])

Notice some features of the betas array:

• There is one row per column in the design X.

• There is one column per column in the data Y.

• The first row of betas contains the parameters corresponding to the first column of the design. The first row therefore contains the slope parameters.

• The second row contains the parameters corresponding to the second column of the design. The second row therefore contains the intercept parameters.

• The first column of betas contains the slope, intercept for the first data vector y_0.

Of course this trick will work for any number of columns of Y.

This trick is important because it allows us to estimate the same design on huge number of data vectors very efficiently. In imaging, this is useful because we can arrange our image data in a one row per volume, one column per voxel array, and use technique here for very fast estimation.