2.2. Inverse Problem Using Restricted Maximum Likelihood Method

The inverse problems (EEG only, DOT only, and joint, detailed in Sec. 2.3.2) were all solved using an ReML algorithm.²⁹ If the forward problem is defined in the linear form,²⁹

where Y is the measurement vector, H is the forward matrix, β is the neuronal source vector, and ϵ is the sensor noise, the algorithm attempts to solve the following optimization problem:

where CN denotes the covariance matrix of measurement noise, CP denotes the covariance matrix of the prior distribution of the neuronal sources, and for some arbitrary matrices X and M, the notation ‖X‖M denotes the weighted norm: ‖X‖M2=XTMX. Such formulation is based on the assumption that both the neuronal sources and the sensor noise follow zero-mean normal distributions, and is derived from the maximum a posteriori estimation. It is worth noting that when CN and CP are both diagonal matrices with equal diagonal elements, the problem reduces to the commonly used Tikhonov regularization.²⁹ In EEG, β represents the electrical activities, Y indicates the scalp EEG recordings, and H is the leadfield matrix calculated using Fieldtrip.²⁵ In DOT, β represents ΔHbO and ΔHb, Y represents ΔOD, and H is the Jacobian calculated using NIRFAST.^4,29

When using ReML, instead of solving directly for CN and CP, structural assumptions can be made on the covariance matrices by rewriting them in forms of linear decomposition, i.e.,

where QN,i and QP,i are the symmetric matrices representing the components to construct the covariance matrices, and ΛN,i and ΛP,i are the coefficients to be estimated from the data. Such decomposition provides one with greater flexibility when making assumptions on the covariances, e.g. different wavelengths in DOT may have different measurement noises, source voxels in different brain regions may have different levels of activities.