Solving the inverse problem

The talon package, provides a way to solve the following optimization problem

\[x^* = \arg\!\min_x \frac{1}{2}\|Ax - y\|_2^2 + \Omega(x)\]

where \(x\) is a vector in \(\mathbb{R}^n\), \(A\) is a linear operator from \(\mathbb{R}^n \to \mathbb{R}^m\) and \(y\) is a vector in \(\mathbb{R}^m\). The functional \(\Omega: \mathbb{R}^n \to \mathbb{R}\) acts as regularization term and must be convex and lower semi-continuous.

The first term of the target functional is devoted to the fitting of the data vector by means of the forward linear operator \(A\) and the coefficient \(x_j\) associated to each atom of \(A\).

Defining regularization term

The possible choices for the regularization term are the following.

• Least Squares

• Non Negativity Constraint

• Structured Sparsity

• Structured Sparsity with Non Negativity

Each of these regularization terms can be defined in talon by calling the talon.regularization function.

Least Squares

Whenever \(\Omega(x) = 0\) for all the admissible values of \(x\), the problem reduces to the classical Least Squares formulation. This is the default regularization term in talon, hence one just needs to call the talon.regularization function as follows.

regterm = talon.regularization()

See an example of this problem in Solve the Least Squares problem.

Non Negativity Constraint

To solve the Non Negative Least Squares (NNLS) problem the regularization term must be the indicator function (in the sense of convex analysis) of the first orthant, namely

\[\Omega(x) = \iota_{\ge 0}(x)\]

which is the function that takes value \(\infty\) whenever \(x\) does not belong to the first orthant. The talon way to obtain such a regularization term is the following.

regterm = talon.regularization(non_negativity=True)

See an example of this problem in Solve the Non Negative Least Squares (NNLS) problem.

Structured Sparsity

To promote sparse solutions, define the group sparsity regularization term

\[\Omega(x) = \lambda \sum_{g\in G} w_g \|x_g\|_2\]

where \(\lambda\) is the regularization parameter, \(w_g\) is the weight associated to each group \(g\), \(x_g\) is the subset of entries of \(x\) corresponding to group \(g\) and \(G\) is the list of groups. See [2011j] for a discussion on the mathematical definition of these groups.

The groups \(g\in G\) must be defined as a list of lists, where each element encodes the indices that define a single group. The weights \(w_g\) associated to each group must be contained in a single numpy array of the same length as \(G\). The following code defines three groups and some standard weight for each of them.

groups = [[0, 2, 5], [1, 3, 4, 6], [7, 8, 9]]
weights = np.array([1.0 / len(g) for g in groups])

Ones the groups, the weights and the regularization parameter are defined, the regularization term can be initialized as follows.

print('Regularization parameter: {}'.format(the_lambda))
print('Number of groups: {}'.format(len(groups)))
print('Number of weights: {}'.format(len(weights)))

regterm = talon.regularization(regularization_parameter=the_lambda,
                               groups=groups, weights=weights)

See an example of this problem at Solve the Group Sparsity problem.

Notice that the standard \(\ell_1\) regularization is a particular case of structure sparsity where there is only one group containing all the admissible indices. Assuming that these indices are \(0\dots n\), the following line of code defines the problem for classical \(\ell_1\) regularization.

groups = [list(range(n))]

See an example of this problem at Solve the Lasso problem and Solve the Non Negative Lasso problem.

Structured Sparsity with Non Negativity

To add the Non Negativity constraint to the Structured Sparsity regularization we just need to set the non_negativity flag as True during the initialization of the regularization term.

regterm = talon.regularization(regularization_parameter=the_lambda,
                               groups=groups, weights=weights,
                               non_negativity=True) # here it is

See an example of this problem at Solve the Non Negative Group Sparsity problem.

Computing the solution

The function devoted to the computation of the solution of the inverse problem is the talon.solve function. It can be called as follows.

linear_operator = # build linear operator
data = # define the data to fit
reg_term = # initialize the regularization term as above

solution = talon.solve(linear_operator=linear_operator,
                       data=data,
                       reg_term=regterm)

The optimization problem is solved with the FISTA+BT algorithm proposed by Beck and Teboulle in [2009b].

See the API documentation for the description of the supplementary optional parameters.

The talon.solve function is a wrapper of the pyunlocbox.solvers.solve function.

Reading the result

The result of the optimization problem is given as a scipy.optimize.OptimizeResult object, which is a dictionary with the following fields.

• x: estimated solution.

• status: attribute of talon.solve.ExitStatus enumeration. If

  \(\text{status} < 1\), the algorithm didn’t converge properly.

• message: string explaining reason for termination.

• fun: value of the objective function at the minimizer.

• nit: number of performed iterations

• reg_param: value of the regularization parameter, if employed.

Examples

Build the ground truth tractogram with two bundles of fibers.

import matplotlib.pyplot as plt
import numpy as np
import talon

from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import interp1d

# Set seed for reproducibility
np.random.seed(1992)

# The number of voxels in each dimension of the output image.
image_size = 25
center = image_size // 2

n_points = int(image_size / 0.01)
t = np.linspace(0, 1, n_points)

# Generate the ground truth tractogram.
tractogram = []
n_streamlines_per_bundle = 50

horizontal_points = np.array([[0, center, center],
                             [image_size - 1, center, center]])
horizontal_streamline = interp1d([0, 1], horizontal_points, axis=0)(t)

for k in range(n_streamlines_per_bundle):
    new_streamline = horizontal_streamline.copy()
    new_streamline[:,1] += (np.random.rand(1) - 0.5)
    tractogram.append(new_streamline)

vertical_points = np.array([[center, 0, center],
                           [center, image_size - 1, center]])
vertical_streamline = interp1d([0, 1], vertical_points, axis=0)(t)

for k in range(n_streamlines_per_bundle):
    new_streamline = vertical_streamline.copy()
    new_streamline[:,0] += (np.random.rand(1) - 0.5)
    tractogram.append(new_streamline)

Show the ground truth tractogram.

fig = plt.figure(figsize=(5, 5), dpi=150)
ax = fig.add_subplot(111, projection='3d')

for streamline in tractogram:
    ax.plot(streamline[:,0], streamline[:,1], streamline[:,2], 'r',
            linewidth=0.1)

ax.plot(horizontal_streamline[:,0],
        horizontal_streamline[:,1],
        horizontal_streamline[:,2], 'k')
ax.plot(vertical_streamline[:,0],
        vertical_streamline[:,1],
        vertical_streamline[:,2], 'k')
ax.view_init(90,90)
ax.set_zticks([])
plt.title('Ground truth tractogram')
plt.show()

You should see the following image:

Generate the corresponding linear operator and the streamline density.

directions = talon.utils.directions(1000)
generators = np.ones((len(directions), 1))
image_shape = (image_size,) * 3
indices, lengths = talon.voxelize(tractogram, directions, image_shape)
linear_operator = talon.operator(generators, indices, lengths)

data = linear_operator @ np.ones(linear_operator.shape[1], dtype=np.float64)
image = data.reshape(image_shape)

Plot the density of the ground truth streamlines

plt.figure(figsize=(5, 5), dpi=150)
plt.imshow(image[:, :, center])
plt.colorbar(shrink=0.8)
plt.title('Ground truth density of streamlines')
plt.show()

You should see the following image:

Add a diagonal bundle of false positives.

diagonal_points = np.array([[0, center, center],
                           [center, image_size - 1, center]])
diagonal_streamline = interp1d([0, 1], diagonal_points, axis=0)(t)

for k in range(n_streamlines_per_bundle):
    new_streamline = diagonal_streamline.copy()
    new_streamline[:,0] += (np.random.rand(1) - 0.5)
    new_streamline[:,1] += (np.random.rand(1) - 0.5)
    tractogram.append(new_streamline)

Visualize the new tractogram.

fig = plt.figure(figsize=(5, 5), dpi=150)
ax = fig.add_subplot(111, projection='3d')

for streamline in tractogram:
    ax.plot(streamline[:,0], streamline[:,1], streamline[:,2], 'r', linewidth=0.1)

ax.plot(horizontal_streamline[:,0],
        horizontal_streamline[:,1],
        horizontal_streamline[:,2], 'k')
ax.plot(vertical_streamline[:,0],
        vertical_streamline[:,1],
        vertical_streamline[:,2], 'k')
ax.plot(diagonal_streamline[:,0],
        diagonal_streamline[:,1],
        diagonal_streamline[:,2], 'k')

ax.view_init(90,90)
ax.set_zticks([])
plt.title('Tractogram with supplementary bundle')
plt.show()

You should see the following image:

Define the linear operator of the tractogram.

indices, lengths = talon.voxelize(tractogram, directions, image_shape)
linear_operator = talon.operator(generators, indices, lengths)

Solve the Least Squares problem

solution = talon.solve(linear_operator=linear_operator, data=data,
                       verbose='NONE')

print('\nLeast Squares solution')
print('Success: {}'.format(solution['success']))
print('Status: {}'.format(solution['status']))
print('Exit criterion: {}'.format(solution['message']))
print('Number of iterations: {}'.format(solution['nit']))

x = solution['x']
print('Average coefficient of horizontal streamlines: {}'.format(
      np.sum(x[0:n_streamlines_per_bundle])/n_streamlines_per_bundle))
print('Average coefficient of vertical streamlines: {}'.format(
      np.sum(x[n_streamlines_per_bundle:2*n_streamlines_per_bundle])/
      n_streamlines_per_bundle))
print('Average coefficient of diagonal streamlines : {}'.format(
      np.sum(x[2*n_streamlines_per_bundle:3*n_streamlines_per_bundle])/
      n_streamlines_per_bundle))
print('Value at minimizer: {}\n'.format(sum(solution['fun'])))

The output should be the following.

Least Squares solution
Success: True
Status: ExitStatus.ABSOLUTE_TOLERANCE_X
Exit criterion: XTOL
Number of iterations: 145
Average coefficient of horizontal streamlines: 0.9999996764340565
Average coefficient of vertical streamlines: 0.9999996573175529
Average coefficient of diagonal streamlines : 4.908558143242968e-06
Value at minimizer: 7.0157355592255e-07

Solve the Non Negative Least Squares (NNLS) problem

reg_term = talon.regularization(non_negativity=True)
solution = talon.solve(linear_operator=linear_operator, data=data,
                       reg_term=reg_term, verbose='NONE')

print('\nNNLS solution')
print('Success: {}'.format(solution['success']))
print('Status: {}'.format(solution['status']))
print('Exit criterion: {}'.format(solution['message']))
print('Number of iterations: {}'.format(solution['nit']))

x = solution['x']
print('Average coefficient of horizontal streamlines: {}'.format(
      np.sum(x[0:n_streamlines_per_bundle])/n_streamlines_per_bundle))
print('Average coefficient of vertical streamlines: {}'.format(
      np.sum(x[n_streamlines_per_bundle:2*n_streamlines_per_bundle])/
      n_streamlines_per_bundle))
print('Average coefficient of diagonal streamlines : {}'.format(
      np.sum(x[2*n_streamlines_per_bundle:3*n_streamlines_per_bundle])/
      n_streamlines_per_bundle))
print('Value at minimizer: {}\n'.format(sum(solution['fun'])))

The output should be the following.

NNLS solution
Success: True
Status: ExitStatus.ABSOLUTE_TOLERANCE_X
Exit criterion: XTOL
Number of iterations: 25
Average coefficient of horizontal streamlines: 0.9999991567472424
Average coefficient of vertical streamlines: 0.9999991568721199
Average coefficient of diagonal streamlines : 5.0072499918376545e-06
Value at minimizer: 3.620593044727195e-07

Solve the Lasso problem

regpar = 1.0 # regularization parameter a.k.a. the lambda in the formula
groups = []
groups.append([k for k in range(0, len(tractogram))])

weights = np.array([1.0 / np.sqrt(len(g)) for g in groups])

reg_term = talon.regularization(groups=groups, weights=weights,
                                regularization_parameter=regpar)

solution = talon.solve(linear_operator=linear_operator, data=data,
                       reg_term=reg_term, verbose='NONE')
print('\nLasso solution')
print('Success: {}'.format(solution['success']))
print('Status: {}'.format(solution['status']))
print('Exit criterion: {}'.format(solution['message']))
print('Number of iterations: {}'.format(solution['nit']))

x = solution['x']
print('Average coefficient of horizontal streamlines: {}'.format(
      np.sum(x[0: n_streamlines_per_bundle])/n_streamlines_per_bundle))
print('Average coefficient of vertical streamlines: {}'.format(
      np.sum(x[n_streamlines_per_bundle:2*n_streamlines_per_bundle]) /
      n_streamlines_per_bundle))
print('Average coefficient of diagonal streamlines : {}'.format(
      np.sum(x[2 * n_streamlines_per_bundle: 3 * n_streamlines_per_bundle]) /
      n_streamlines_per_bundle))
print('Value at minimizer: {}\n'.format(sum(solution['fun'])))

The output should be the following:

Lasso solution
Success: True
Status: ExitStatus.RELATIVE_TOLERANCE_COST
Exit criterion: RTOL
Number of iterations: 93
Average coefficient of horizontal streamlines: 0.9999926298816814
Average coefficient of vertical streamlines: 0.9999925070704963
Average coefficient of diagonal streamlines : -2.1995490196016877e-05
Value at minimizer: 0.8165122997013363

Solve the Non Negative Lasso problem

reg_term = talon.regularization(non_negativity=True,
                                groups=groups, weights=weights,
                                regularization_parameter=regpar)

solution = talon.solve(linear_operator=linear_operator, data=data,
                       reg_term=reg_term, verbose='NONE')
print('\nNon Negative Lasso solution')
print('Success: {}'.format(solution['success']))
print('Status: {}'.format(solution['status']))
print('Exit criterion: {}'.format(solution['message']))
print('Number of iterations: {}'.format(solution['nit']))

x = solution['x']
print('Average coefficient of horizontal streamlines: {}'.format(
      np.sum(x[0: n_streamlines_per_bundle])/n_streamlines_per_bundle))
print('Average coefficient of vertical streamlines: {}'.format(
      np.sum(x[n_streamlines_per_bundle:2*n_streamlines_per_bundle]) /
      n_streamlines_per_bundle))
print('Average coefficient of diagonal streamlines : {}'.format(
      np.sum(x[2 * n_streamlines_per_bundle: 3 * n_streamlines_per_bundle]) /
      n_streamlines_per_bundle))
print('Value at minimizer: {}\n'.format(sum(solution['fun'])))

The output should be the following:

Non Negative Lasso solution
Success: True
Status: ExitStatus.RELATIVE_TOLERANCE_COST
Exit criterion: RTOL
Number of iterations: 23
Average coefficient of horizontal streamlines: 0.9999914147578718
Average coefficient of vertical streamlines: 0.9999914603196133
Average coefficient of diagonal streamlines : 4.482209580050452e-06
Value at minimizer: 0.8164938196507543

Solve the Group Sparsity problem

groups = []
groups.append([k for k in range(0, n_streamlines_per_bundle)]) # horizontal
groups.append([k for k in range(n_streamlines_per_bundle,
              2 * n_streamlines_per_bundle)]) # vertical
groups.append([k for k in range(2 * n_streamlines_per_bundle,
              3 * n_streamlines_per_bundle)]) # diagonal

weights = np.array([1.0 / np.sqrt(len(g)) for g in groups])

reg_term = talon.regularization(groups=groups, weights=weights,
                                regularization_parameter=regpar)

solution = talon.solve(linear_operator=linear_operator, data=data,
                       reg_term=reg_term, verbose='NONE')
print('\nGroup Sparsity solution')
print('Success: {}'.format(solution['success']))
print('Status: {}'.format(solution['status']))
print('Exit criterion: {}'.format(solution['message']))
print('Number of iterations: {}'.format(solution['nit']))

x = solution['x']
print('Average coefficient of horizontal streamlines: {}'.format(
      np.sum(x[0: n_streamlines_per_bundle])/n_streamlines_per_bundle))
print('Average coefficient of vertical streamlines: {}'.format(
      np.sum(x[n_streamlines_per_bundle:2*n_streamlines_per_bundle]) /
      n_streamlines_per_bundle))
print('Average coefficient of diagonal streamlines : {}'.format(
      np.sum(x[2 * n_streamlines_per_bundle: 3 * n_streamlines_per_bundle]) /
      n_streamlines_per_bundle))
print('Value at minimizer: {}\n'.format(sum(solution['fun'])))

The output should be the following:

Group Sparsity solution
Success: True
Status: ExitStatus.RELATIVE_TOLERANCE_COST
Exit criterion: RTOL
Number of iterations: 64
Average coefficient of horizontal streamlines: 0.9999821712768615
Average coefficient of vertical streamlines: 0.9999823618643954
Average coefficient of diagonal streamlines : 2.2318881330827924e-05
Value at minimizer: 2.000096258909371

Solve the Non Negative Group Sparsity problem

reg_term = talon.regularization(groups=groups, weights=weights,
                                non_negativity=True,
                                regularization_parameter=regpar)

solution = talon.solve(linear_operator=linear_operator, data=data,
                       reg_term=reg_term, verbose='NONE')
print('\nNon Negative Group Sparsity solution')
print('Success: {}'.format(solution['success']))
print('Status: {}'.format(solution['status']))
print('Exit criterion: {}'.format(solution['message']))
print('Number of iterations: {}'.format(solution['nit']))

x = solution['x']
print('Average coefficient of horizontal streamlines: {}'.format(
      np.sum(x[0: n_streamlines_per_bundle])/n_streamlines_per_bundle))
print('Average coefficient of vertical streamlines: {}'.format(
      np.sum(x[n_streamlines_per_bundle:2*n_streamlines_per_bundle]) /
      n_streamlines_per_bundle))
print('Average coefficient of diagonal streamlines : {}'.format(
      np.sum(x[2 * n_streamlines_per_bundle: 3 * n_streamlines_per_bundle]) /
      n_streamlines_per_bundle))
print('Value at minimizer: {}\n'.format(sum(solution['fun'])))

The output should be the following:

Non Negative Group Sparsity solution
Success: True
Status: ExitStatus.RELATIVE_TOLERANCE_COST
Exit criterion: RTOL
Number of iterations: 22
Average coefficient of horizontal streamlines: 0.9999825264666186
Average coefficient of vertical streamlines: 0.9999825878147537
Average coefficient of diagonal streamlines : 0.0
Value at minimizer: 1.9999822314331122

References

2009b

Beck, Amir, and Marc Teboulle. “A fast iterative shrinkage-thresholding algorithm for linear inverse problems.” SIAM journal on imaging sciences 2.1 (2009): 183-202.

2011j

Jenatton, Rodolphe, et al. “Proximal methods for hierarchical sparse coding.” Journal of Machine Learning Research 12.Jul (2011): 2297-2334.