Adapting source grid parameters to improve the condition of the magnetostatic linear inverse problem of estimating nanoparticle distributions

Abstract

The problem of estimating magnetic nanoparticle distributions from magnetorelaxometric measurements is addressed here. The objective of this work was to identify source grid parameters that provide a good condition for the related linear inverse problem. The parameters investigated here were the number of sources, the extension of the source grid, and the source direction. A new measure of the condition, the ratio between the largest and mean singular value of the lead field matrix, is proposed. Our results indicated that the source grids should be larger than the sensor area. The sources and, consequently, the magnetic excitation field, should be directed toward the Z-direction. For underdetermined linear inverse problems, such as in our application, the number of sources affects the condition to a relatively small degree. Overdetermined magnetostatic linear inverse problems, however, benefit from a reduction in the number of sources, which considerably improves the condition. The adapted source grids proposed here were used to estimate the magnetostatic dipole from simulated data; the L2-norm, residual, and distances between the estimated and simulated sources were significantly reduced.

Appendix

1.1 Accuracy of the condition number κ

The singular value decomposition of the lead field matrices in this application was accomplished using the LAPACK’s DGESVD routine. Anderson et al. (section 4.9.1 of [1]) show that the DGESVD produces deviations in the computed SVs of $\hat{\sigma}_i$ relative to the true values $\sigma_i,$ bounded by

$$ || \hat{\sigma}_i - \sigma_i || \leq \epsilon \sigma_1 . $$

(8)

The floating-point relative accuracy $\epsilon$ was 2.2e−16 on the 64-bit systems used for our computations. Therefore, high relative accuracies of $\hat{\sigma}_i$ are only obtained for SVs close to $\sigma_1.$ Because the smallest SV $\sigma_n$ is very small in relation to $\sigma_1$ in this and many other applications, even low absolute errors in $\sigma_n$ indicate high relative errors and lead to high absolute errors in $\sigma_n^{-1}$ and $\kappa.$

Considering the definition (2), the accuracy of $\kappa$ crucially depends on $\sigma_n.$ Using the inequality (8), we obtain, for the relative error of $\sigma_n$

$$\frac{|| \hat{\sigma}_n - \sigma_n ||}{\sigma_n} \leq \frac{\epsilon \sigma_1}{\sigma_n} = \epsilon \kappa.$$

(9)

Consequently, the values of $\kappa$ may be inaccurate, even in the order of magnitude, for $\kappa \geq \epsilon^{-1}.$ Demmel [8] and Higham [15] stated that the condition number for computing the condition number is the condition number.

For special classes of matrices, it is feasible to compute all SVs, including the tiny SVs, with a high relative accuracy. An overview of this topic is provided by Demmel et al. [9]. This is not possible for general dense matrices, such as the lead field operators employed in our application.

1.2 Condition of the TSVD-regularized linear inverse problem

To facilitate stable linear inverse solutions, regularization methods are applied to improve the condition. For example, when computing a linear inverse solution using the TSVD approach (7) and a regularization parameter r, with $1\leq r \leq n,$ all singular components with SVs smaller than $\sigma_r$ are omitted. Depending on the parameter r, the TSVD regularization changes the linear IP and causes a regularization error, which should be smaller than the error that results from an inferior condition.

To measure the condition of the TSVD-regularized linear inverse problem for a lead field matrix L, the definitions of $\kappa$ and ρ, (2) and (3), are modified to consider the largest r SVs only:

$$ \kappa_{\rm tsvd}({\mathbf {L}},r) = \sigma_1({\mathbf {L}}) \sigma_r({\mathbf {L}})^{-1} , $$

(10)

$$ \rho_{\rm tsvd}({\mathbf {L}},r) = \frac{\sigma_1({\mathbf {L}})} {1/r \sum_{i=1}^r \sigma_i({\mathbf {L}})} = \left( \frac{1}{r} \sum_{i=1}^r \frac{\sigma_i({\mathbf {L}})} {\sigma_1({\mathbf {L}})} \right)^{-1} . $$

(11)

From (2) and (10), we obtain for any r between 1 and n

$$ 1 \leq \kappa_{\rm tsvd}({\mathbf {L}},r) \leq \kappa({\mathbf {L}}) , $$

(12)

and, from (3) and (11),

$$ 1 \leq \rho_{\rm tsvd}({\mathbf{L}},r) \leq \rho({\mathbf {L}}) . $$

(13)

The property (12) follows because $\sigma_1\geq\sigma_r\geq\sigma_n.$ A proof for (13) is shown in Sect. “Proof of $1\leq \rho_{\rm tsvd} \leq \rho$” of the Appendix. According to (12) and (13), $\kappa$ and ρ are the least upper bounds for $\kappa_{\rm tsvd}$ and $\rho_{\rm tsvd}.$ In practice, the optimal choice of r for a matrix L and the resulting values of $\kappa_{\rm tsvd}({\mathbf{{L}}},r)$ and $\rho_{\rm tsvd}({\mathbf{{L}}},r)$ can vary widely.

1.3 Proof of $1\leq \rho_{\rm tsvd} \leq \rho$

To prove $1 \leq \rho_{\rm tsvd}({\mathbf {L}},r) \leq \rho({\mathbf {L}}),$ we first show that $\rho_{\rm tsvd}({\mathbf {L}},r)$ is smaller or equal $\rho({\mathbf {L}})$ for all regularization parameters r between 1 and n:

$$\begin{aligned} &\left( \frac{1}{r} \sum_{i=1}^r \frac{\sigma_i({\mathbf {L}})} {\sigma_1({\mathbf {L}})} \right)^{-1} \leq \left( \frac{1}{n} \sum_{i=1}^n \frac{\sigma_i({\mathbf {L}})} {\sigma_1({\mathbf {L}})} \right)^{-1}\\ &\quad\Longleftrightarrow \frac{1}{r} \sum_{i=1}^r \frac{\sigma_i({\mathbf {L}})} {\sigma_1({\mathbf {L}})} \geq \frac{1}{n} \sum_{i=1}^n \frac{\sigma_i({\mathbf {L}})} {\sigma_1({\mathbf {L}})}\\&\quad \Longleftrightarrow n \sum_{i=1}^r \sigma_i({\mathbf {L}}) \geq r \sum_{i=1}^n \sigma_i({\mathbf {L}})\\ &\quad\Longleftrightarrow r \sum_{i=1}^r \sigma_i({\mathbf {L}}) +(n-r) \underbrace{\sum_{i=1}^r \sigma_i({\mathbf {L}})}_{\geq r\times\sigma_r({\mathbf {L}})} \ \geq r \sum_{i=1}^r \sigma_i({\mathbf {L}}) + r \underbrace{\sum_{i=r+1}^n \sigma_i({\mathbf {L}})}_{\leq (n-r)\times\sigma_{r+1}({\mathbf {L}})}\\&\quad \Longleftarrow \, (n-r)\times r\times\sigma_r({\mathbf {L}}) \geq r\times(n-r)\times\sigma_{r+1}({\mathbf {L}})\\&\quad \Longleftrightarrow \sigma_r({\mathbf {L}}) \geq \sigma_{r+1}({\mathbf {L}}) \end{aligned} $$

which is provided by the singular value decomposition of L. Second, from $\sigma_{1}({\mathbf {L}})\geq\sigma_{i}({\mathbf {L}})$ follows $r \times \sigma_1({\mathbf {L}}) \geq \sum_{i=1}^r \sigma_i({\mathbf {L}})$ and $1 \leq \rho_{\rm tsvd}({\mathbf {L}},r).$ $\square$