Coverage for local_installation_linux/mumott/methods/utilities/preconditioning.py: 100%
41 statements
« prev ^ index » next coverage.py v7.3.2, created at 2025-05-05 21:21 +0000
1import logging
3import numpy as np
4from numpy.typing import NDArray
6from mumott.methods.projectors.base_projector import Projector
7from mumott.methods.basis_sets.base_basis_set import BasisSet
9logger = logging.getLogger(__name__)
12def get_largest_eigenvalue(basis_set: BasisSet,
13 projector: Projector,
14 weights: NDArray[float] = None,
15 preconditioner: NDArray[float] = None,
16 niter: int = 5,
17 seed: int = None) -> float:
18 r"""
19 Calculate the largest eigenvalue of the matrix :math:`A^{T}*A` using the power method.
20 Here, :math:`A` is the linear forward model defined be the input
21 :class:`Projector <mumott.method.projectors.base_projector.Projector>` and
22 :class:`BasisSet <mumott.method.bais_sets.base_basis_set.BasisSet>`.
23 The largest eigenvalue can be used to set a safe step size for various optimizers.
25 Parameters
26 ----------
27 basis_set
28 basis set object
29 projector
30 projector object
31 weights
32 Weights which will be applied to the residual.
33 Default is ``None``.
34 preconditioner
35 A preconditioner which will be applied to the gradient.
36 Default is ``None``.
37 niter
38 number of iterations. Default is 5
39 seed
40 Seed for random generation as starting state. Used for testing.
42 Returns.
43 -------
44 An estimate of the matrix norm (largest singular value)
45 """
46 shape = (*projector.volume_shape, len(basis_set))
48 if seed is not None:
49 np.random.seed(seed)
50 if weights is None:
51 weights = 1
52 if preconditioner is None:
53 preconditioner = 1
55 x = np.random.normal(loc=0.0, scale=1.0, size=shape).astype(projector.dtype)
57 for _ in range(niter):
58 x = x / np.sqrt(np.sum(x**2))
59 y = basis_set.forward(projector.forward(x)) * weights
60 x = projector.adjoint(basis_set.gradient(y).astype(projector.dtype)) * \
61 preconditioner
63 return np.sqrt(np.sum(x**2))
66def get_tensor_sirt_preconditioner(projector: Projector,
67 basis_set: BasisSet,
68 cutoff: float = 0.1) -> NDArray[float]:
69 r""" Retrieves the :term:`SIRT` `preconditioner <https://en.wikipedia.org/wiki/Preconditioner>`_
70 for tensor representations, which can be used together
71 with the tensor :term:`SIRT` weights to condition the
72 gradient of tensor tomographic calculations for faster convergence and scaling of the
73 step size for gradient descent.
75 Notes
76 -----
77 The preconditioner is computed similarly as in the scalar case, except the underlying
78 system of equations is
80 .. math::
81 P_{ij} X_{jk} U_{kl} = Y_{il}
83 where :math:`X_{jk}` is a vector of tensor-valued voxels, and :math:`Y_{il}`
84 is the projection into measurement space. The preconditioner then corresponds to
86 .. math::
87 C_{jjkk} = \frac{I(n_j) \otimes I(n_k)}{\sum_{i,l} P_{ij} U_{kl}},
89 where :math:`I(n_j)` is the identity matrix of the same size as :math:`X_j`. Similarly,
90 the weights (of :func:`~.get_tensor_sirt_weights`) are computed as
92 .. math::
93 W_{iill} = \frac{I(n_j) \otimes I(n_k)}{\sum_{j, k} P_{ij} U_{kl}}.
95 Here, any singularities in the system (e.g., where :math:`\sum_j P_{ij} = 0`) can be masked out
96 by setting the corresponding weight to zero.
97 We thus end up with a weighted least-squares system
99 .. math::
100 \text{argmin}_X(\Vert W_{iill}(P_{ij} X_{jk} U_{kl} - D_{il})\Vert^2_2),
102 where :math:`D_{il}` is some data.
103 This problem can be solved iteratively by preconditioned gradient descent,
105 .. math::
106 X_j^{k + 1} = X_j^k - C_{jjkk}P_{ji}^TW_{iill}(P_ij X_{jk}^k U_{kl} - D_{il}).
108 Parameters
109 ----------
110 projector
111 A :ref:`projector <projectors>` object which is used to calculate the weights.
112 The computation of the weights is based on the geometry attached to the projector.
113 basis_set
114 A :ref:`basis set <basis_sets>` object which is used to calculate the weights.
115 The tensor-space-to-detector-space transform uses the basis set.
116 Should use a local representation.
117 cutoff
118 The minimal number of rays that need to map to a voxel for it
119 to be considered valid. Default is ``0.1``. Invalid voxels will
120 be masked from the preconditioner.
121 """
122 sirt_projections = np.ones((projector.number_of_projections,) +
123 projector.projection_shape +
124 (basis_set.probed_coordinates.vector.shape[1],),
125 dtype=projector.dtype)
126 inverse_preconditioner = projector.adjoint(basis_set.gradient(sirt_projections).astype(projector.dtype))
127 mask = np.ceil(inverse_preconditioner > projector.dtype(cutoff)).astype(projector.dtype)
128 return mask * np.reciprocal(np.fmax(inverse_preconditioner, cutoff))
131def get_sirt_preconditioner(projector: Projector, cutoff: float = 0.1) -> NDArray[float]:
132 r""" Retrieves the :term:`SIRT` `preconditioner <https://en.wikipedia.org/wiki/Preconditioner>`_,
133 which can be used together
134 with the :term:`SIRT` weights to condition the
135 gradient of tomographic calculations for faster convergence and scaling of the
136 step size for gradient descent.
138 Notes
139 -----
140 The preconditioner normalizes the gradient according to the number
141 of data points that map to each voxel in the computation of
142 the projection adjoint. This preconditioner scales and conditions
143 the gradient for better convergence. It is best combined with the :term:`SIRT`
144 weights, which normalize the residual for the number of voxels.
145 When used together, they condition the reconstruction sufficiently
146 well that a gradient descent optimizer with step size unity can arrive
147 at a good solution. Other gradient-based solvers can also benefit from this
148 preconditioning.
150 In addition, the calculation of these preconditioning terms makes it easy to identify
151 regions of the volume or projections that are rarely probed, allowing them to be
152 masked from the solution altogether.
154 If the projection operation is written in sparse matrix form as
156 .. math::
157 P_{ij} X_{j} = Y_i
159 where :math:`P_{ij}` is the projection matrix, :math:`X_j` is a vector of voxels, and :math:`Y_i`
160 is the projection, then the preconditioner can be understood as
162 .. math::
163 C_{jj} = \frac{I(n_j)}{\sum_i P_{ij}}
165 where :math:`I(n_j)` is the identity matrix of the same size as :math:`X_j`. Similarly,
166 the weights (of :func:`~.get_sirt_weights`) are computed as
168 .. math::
169 W_{ii} = \frac{I(n_i)}{\sum_j P_{ij}}.
171 Here, any singularities in the system (e.g., where :math:`\sum_j P_{ij} = 0`) can be masked out
172 by setting the corresponding weight to zero.
173 We thus end up with a weighted least-squares system
175 .. math::
176 \text{argmin}_X(\Vert W_{ii}(P_{ij}X_{j} - D_{i})\Vert^2_2)
178 where :math:`D_{i}` is some data, which we can solve iteratively by preconditioned gradient descent,
180 .. math::
181 X_j^{k + 1} = X_j^k - C_{jj}P_{ji}^TW_{ii}(P_ij X_j^k - D_i)
183 As mentioned, we can add additional regularization terms, and because the preconditioning
184 scales the problem appropriately, computing an optimal step size is not a requirement,
185 although it can speed up the solution. This establishes a very flexible system, where
186 quasi-Newton solvers such as :term:`LBFGS` can be seamlessly combined with less restrictive
187 gradient descent methods.
189 A good discussion of the algorithmic properties of :term:`SIRT` can be found in
190 `this article by Gregor et al. <https://doi.org/10.1109%2FTCI.2015.2442511>`_,
191 while `this article by van der Sluis et al. <https://doi.org/10.1016/0024-3795(90)90215-X>`_
192 discusses :term:`SIRT` as a least-squares solver in comparison to the
193 conjugate gradient (:term:`CG`) method.
195 Parameters
196 ----------
197 projector
198 A :ref:`projector <projectors>` object which is used to calculate the weights.
199 The computation of the weights is based on the geometry attached to the projector.
200 cutoff
201 The minimal number of rays that need to map to a voxel for it
202 to be considered valid. Default is ``0.1``. Invalid voxels will
203 be masked from the preconditioner.
204 """
205 sirt_projections = np.ones((projector.number_of_projections,) +
206 projector.projection_shape +
207 (1,), dtype=projector.dtype)
208 inverse_preconditioner = projector.adjoint(sirt_projections)
209 mask = np.ceil(inverse_preconditioner > projector.dtype(cutoff)).astype(projector.dtype)
210 return mask * np.reciprocal(np.fmax(inverse_preconditioner, cutoff))
213def get_tensor_sirt_weights(projector: Projector, basis_set: BasisSet, cutoff: float = 0.1) -> NDArray[float]:
214 """ Retrieves the tensor :term:`SIRT` weights, which can be used together with the
215 :term:`SIRT` preconditioner to weight the
216 residual of tensor tomographic calculations for faster convergence and
217 scaling of the step size for gradient descent.
219 Notes
220 -----
221 See :func:`~.get_tensor_sirt_preconditioner` for theoretical details.
223 Parameters
224 ----------
225 projector
226 A :ref:`projector <projectors>` object which is used to calculate the weights.
227 The computation of the weights is based on the geometry attached to the projector.
228 basis_set
229 A :ref:`basis set <basis_sets>` object which is used to calculate the weights.
230 The tensor-space-to-detector-space transform uses the basis set.
231 Should use a local representation.
232 cutoff
233 The minimal number of voxels that need to map to a point for it
234 to be considered valid. Default is ``0.1``. Invalid pixels will be
235 masked.
236 """
237 sirt_field = np.ones(projector.volume_shape +
238 (len(basis_set),), dtype=projector.dtype)
239 inverse_weights = projector.forward(sirt_field)
240 inverse_weights = basis_set.forward(inverse_weights)
241 mask = np.ceil(inverse_weights > projector.dtype(cutoff)).astype(projector.dtype)
242 return mask * np.reciprocal(np.fmax(inverse_weights, cutoff))
245def get_sirt_weights(projector: Projector, cutoff: float = 0.1) -> NDArray[float]:
246 """ Retrieves the :term:`SIRT` weights, which can be used together with the
247 :term:`SIRT` preconditioner to weight the
248 residual of tomographic calculations for faster convergence and
249 scaling of the step size for gradient descent.
251 Notes
252 -----
253 See :func:`~.get_sirt_preconditioner` for theoretical details.
255 Parameters
256 ----------
257 projector
258 A :ref:`projector <projectors>` object which is used to calculate the weights.
259 The computation of the weights is based on the geometry attached to the projector.
260 cutoff
261 The minimal number of voxels that need to map to a point for it
262 to be considered valid. Default is ``0.1``. Invalid pixels will be
263 masked.
264 """
265 sirt_field = np.ones(projector.volume_shape +
266 (1,), dtype=projector.dtype)
267 inverse_weights = projector.forward(sirt_field)
268 mask = np.ceil(inverse_weights > projector.dtype(cutoff)).astype(projector.dtype)
269 return mask * np.reciprocal(np.fmax(inverse_weights, cutoff))