pct_tools.recon

Steepest descent on 1/2 ||A x - b||^2 for a system stored as blocks on disk.

The system matrix A and the measurement vector b are split into row blocks (A_i, b_i), one pair of files per projection, written by pct_tools.matrix and the preprocessing. Each call here streams through every block in parallel and returns one quantity for the whole system: one descent step, or the squared residual norm. The iteration itself, including the step size (see pct_tools.lipschitz) and the stopping rule, is left to the caller.

pct_tools.recon.gradient_descent(x, step_size, matrix_filenames, vector_filenames, *, verbosity_level=0, n_threads=1)[source]

Perform one steepest-descent step on 1/2 ||A x - b||^2.

The system is split into blocks (A_i, b_i) that are stored on disk and processed in parallel; the result is x - step_size * sum_i A_i^T (A_i x - b_i).

Parameters:

  • x (NDArray[np.float32]) – The current estimate, a float32 array of shape (n_voxels,) where n_voxels is the number of columns of A.

  • step_size (float) – The step size to be used in the descent, typically 1 / L with L the Lipschitz constant estimated via pct_tools.lipschitz.

  • matrix_filenames (Sequence[str]) – The filenames of the compressed matrix blocks.

  • vector_filenames (Sequence[str]) – The filenames of the matching b-vector blocks.

  • verbosity_level (int) – The verbosity level.

  • n_threads (int) – The number of threads to use.

Returns:

The updated vector after one step of steepest descent, with the same shape as x.

Raises:

TypeError – If x is not a float32 array.

Return type:

NDArray[np.float32]

pct_tools.recon.squared_residuals(x, matrix_filenames, vector_filenames, *, verbosity_level=0, n_threads=1)[source]

Compute the squared residual norm ||A x - b||^2 for a given vector.

Parameters:

  • x (NDArray[np.float32]) – The current estimate, a float32 array of shape (n_voxels,).

  • matrix_filenames (Sequence[str]) – The filenames of the compressed matrix blocks.

  • vector_filenames (Sequence[str]) – The filenames of the matching b-vector blocks.

  • verbosity_level (int) – The verbosity level.

  • n_threads (int) – The number of threads to use.

Returns:

The squared residual norm.

Raises:

TypeError – If x is not a float32 array.

Return type:

float