.. _inverse-problem: =========================== Solving the inverse problem =========================== The ``talon`` package, provides a way to solve the following optimization problem .. math:: x^* = \arg\!\min_x \frac{1}{2}\|Ax - y\|_2^2 + \Omega(x) where :math:`x` is a vector in :math:`\mathbb{R}^n`, :math:`A` is a linear operator from :math:`\mathbb{R}^n \to \mathbb{R}^m` and :math:`y` is a vector in :math:`\mathbb{R}^m`. The functional :math:`\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 :math:`A` and the coefficient :math:`x_j` associated to each atom of :math:`A`. Defining regularization term ---------------------------- The possible choices for the regularization term are the following. * :ref:`least-squares` * :ref:`non-negativity-constraint` * :ref:`structured-sparsity` * :ref:`structured-sparsity-non-negativity` Each of these regularization terms can be defined in `talon` by calling the :code:`talon.regularization` function. .. _least-squares: Least Squares ^^^^^^^^^^^^^ Whenever :math:`\Omega(x) = 0` for all the admissible values of :math:`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 :code:`talon.regularization` function as follows. .. code:: python regterm = talon.regularization() See an example of this problem in :ref:`solve-least-squares`. .. _non-negativity-constraint: 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 .. math:: \Omega(x) = \iota_{\ge 0}(x) which is the function that takes value :math:`\infty` whenever :math:`x` does not belong to the first orthant. The `talon` way to obtain such a regularization term is the following. .. code:: python regterm = talon.regularization(non_negativity=True) See an example of this problem in :ref:`solve_nnls`. .. _structured-sparsity: Structured Sparsity ^^^^^^^^^^^^^^^^^^^ To promote sparse solutions, define the group sparsity regularization term .. math:: \Omega(x) = \lambda \sum_{g\in G} w_g \|x_g\|_2 where :math:`\lambda` is the regularization parameter, :math:`w_g` is the weight associated to each group :math:`g`, :math:`x_g` is the subset of entries of :math:`x` corresponding to group :math:`g` and :math:`G` is the list of groups. See [2011j]_ for a discussion on the mathematical definition of these groups. The groups :math:`g\in G` must be defined as a list of lists, where each element encodes the indices that define a single group. The weights :math:`w_g` associated to each group must be contained in a single numpy array of the same length as :math:`G`. The following code defines three groups and some standard weight for each of them. .. code:: python 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. .. code:: python 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 :ref:`solve-gs`. Notice that the standard :math:`\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 :math:`0\dots n`, the following line of code defines the problem for classical :math:`\ell_1` regularization. .. code:: python groups = [list(range(n))] See an example of this problem at :ref:`solve-lasso` and :ref:`solve-nnlasso`. .. _structured-sparsity-non-negativity: Structured Sparsity with Non Negativity ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ To add the Non Negativity constraint to the Structured Sparsity regularization we just need to set the :code:`non_negativity` flag as :code:`True` during the initialization of the regularization term. .. code:: python regterm = talon.regularization(regularization_parameter=the_lambda, groups=groups, weights=weights, non_negativity=True) # here it is See an example of this problem at :ref:`solve-nngs`. .. _compute-solution: Computing the solution ---------------------- The function devoted to the computation of the solution of the inverse problem is the :code:`talon.solve` function. It can be called as follows. .. code:: python 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 :code:`talon.solve` function is a wrapper of the :code:`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 :math:`\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. .. code:: python 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. .. code:: python 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: .. image:: img/gt_tractogram.png :align: center :width: 50em :height: 50em Generate the corresponding linear operator and the streamline density. .. code:: python 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 .. code:: python 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: .. image:: img/gt_density.png :align: center :width: 50em :height: 50em Add a diagonal bundle of false positives. .. code:: python 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. .. code:: python 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: .. image:: img/tractogram_with_fp.png :align: center :width: 50em :height: 50em Define the linear operator of the tractogram. .. code:: python indices, lengths = talon.voxelize(tractogram, directions, image_shape) linear_operator = talon.operator(generators, indices, lengths) .. _solve-least-squares: Solve the Least Squares problem ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code:: python 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. .. code:: 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_nnls: Solve the Non Negative Least Squares (NNLS) problem ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code:: python 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. .. code:: 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-lasso: Solve the Lasso problem ^^^^^^^^^^^^^^^^^^^^^^^ .. code:: python 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: .. code:: 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-nnlasso: Solve the Non Negative Lasso problem ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code:: python 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: .. code:: 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-gs: Solve the Group Sparsity problem ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code:: python 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: .. code:: 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-nngs: Solve the Non Negative Group Sparsity problem ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. code:: python 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: .. code:: 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.