The linear regression function is in the form
(1)
where is the muscle activity (regressand) to be reconstructed and

is the design matrix with the time series of a set of p voxels from the selected ROI in the 'regression set'. The parameters will be computed to realize the best mapping of the voxels to the muscle activity while making the function sparse. is the noise for n brain scans recorded in an experiment.
To favour a sparse representation, a Laplacian prior is chosen for β such that

where and α is the parameter of this density, then the maximum a posteriori estimate of β is given by
(2)
The term with the Euclidean norm || ||₂ expresses the least-square approximation of the regress and, while the term || ||₁is to produce a sparse representation. The || ||₁norm, corresponding to the absolute value, has a larger gradient than the || ||₂ norm in the vicinity of the minimum and thus should converge fast to zero, i.e. to a sparse representation. This criterion was named as least absolute shrinkage and selection operator or LASSO (Tibshirani 1996).

The Laplacian prior is equivalent to a two level hierarchical Bayes model with zero mean Gaussian priors with independent exponentially distributed variances (Figueiredo 2003; Lange and Sinsheimer 1993). This allows us to use the expectation- maximization (EM) algorithm to implement the LASSO criterion (Osborne et al. 2000, Tibshirani 1996) in two steps:

In our regression (4) and (5) are used to make β sparse, determined by the value of parameter σ. The optimal value of σ was determined by an adaptive process as explained in next section.

To achieve a good generalization of a model or mapping, it is necessary to control the complexity of the learnt function. If the function is too complex then it will reproduce unimportant details of the training set, thus over-fitting the training data. Such a function will fit the training data well but will be unable to predict a different test data well. Conversely, an overly simple function will not be able to capture the true mapping between the regressor and the regressand, thus under-fitting the training data. In our case the number of voxels chosen during training decides the ability of the model to generalize well. The typical behaviour of our data set is represented in (Fig. 1). The regression algorithm has one free parameter (σ) to adjust the number of voxels used during training. We utilize over-fitting, detected by the performances on a representative portion of the test set, the selection set, as a criterion to recursively select an optimal value of σ.

The EM algorithm loop starts with σ =1. The evaluated mapping vector σ is tested on the selection set and the coefficient of determination (R²) is evaluated. A binary search is used to determine the σ with minimal R².

The performance has been quantified using the coefficient of determination , where y_r represents the muscle regressand value, _r its average and y_c the reconstructed value.

• The regressand is the physical parameter that will be reconstructed/ decoded from the brain activity. We reconstructed the muscle EMG and thus used smoothed EMG activity as regressand.
• The regressand is re-sampled at the same frequency as the repetition time (TR) of the fMRI scan.
• This data is used in the training (yt), test (yr) and selection (ys)

• The brain images from the collected in different sessions should be realigned to the first scan of session 1 to ensure there was no displacement among scans from the different sessions.
• For each subject a region of interest (ROI), may be defined for reconstruction
• The voxel time series from the ROI (or entire brain) is extracted and band-pass filtered between 0.003 Hz (300 sec) and 0.2 Hz (5 sec) before the regression analysis. The low frequency cut-off value will remove any drift in the fMRI signal while the high cut-off of 0.2 Hz will reduce noise constrained by the frequency content of haemodynamic function (hrf). This preprocessed data was used as the regressor in the training (Xt), test (Xr) and selection (Xs) analysis.

• Check the value of R2₁ for a particular sigma (say σ₁=1) using step 4 and step 5.
• Change the sigma, (say σ₂=0.1) and calculate R2₂.
• If R2₂ > R2₁, then next choose s₂ and try values around σ₂.
• Check for sigma values, say σ₃=0.5 and σ₄=0.05 to get R2₃ and R2₄.
• If R2₃>R2₄ then choose σ₃ and try values around σ₃ like step IV else try values around σ₄.
• Continue like in steps IV and V till convergence of R₂ is achieved. This corresponds to the optimal σ_opt (sigma_opt) and gives the optimal beta (beta_opt) from step 4.

• Sparse Regression Toolbox

• sp_regress.m
• reconstruct.m
• data_sample.mat
• Sparse_reg.doc