diff --git a/.github/workflows/pytest_ubuntu.yml b/.github/workflows/pytest_ubuntu.yml index ffad0c04f..babfb1d9f 100644 --- a/.github/workflows/pytest_ubuntu.yml +++ b/.github/workflows/pytest_ubuntu.yml @@ -89,7 +89,7 @@ jobs: - name: Run doc tests if: matrix.group == 1 run: | - pytest --ignore-glob="*test_*.py" --doctest-modules scico + pytest --ignore-glob="*test_*.py" --ignore=scico/linop/xray --doctest-modules scico pytest --doctest-glob="*.rst" docs coverage: diff --git a/data b/data index 199593151..cb97dd02c 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 199593151ab6317207a2a63ed471f86c882b92e2 +Subproject commit cb97dd02cfd96659804c74e7a5f9185440cec3ce diff --git a/docs/source/examples.rst b/docs/source/examples.rst index da89c8289..41a3df9e1 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -38,7 +38,6 @@ Computed Tomography examples/ct_astra_odp_train_foam2 examples/ct_astra_unet_train_foam2 examples/ct_projector_comparison - examples/ct_multi_cs_tv_admm examples/ct_multi_tv_admm Deconvolution @@ -137,6 +136,7 @@ Total Variation examples/ct_astra_3d_tv_admm examples/ct_astra_3d_tv_padmm examples/ct_astra_weighted_tv_admm + examples/ct_multi_tv_admm examples/ct_svmbir_tv_multi examples/deconv_circ_tv_admm examples/deconv_tv_admm @@ -208,6 +208,7 @@ ADMM examples/ct_tv_admm examples/ct_astra_3d_tv_admm examples/ct_astra_weighted_tv_admm + examples/ct_multi_tv_admm examples/ct_svmbir_tv_multi examples/ct_svmbir_ppp_bm3d_admm_cg examples/ct_svmbir_ppp_bm3d_admm_prox diff --git a/docs/source/style.rst b/docs/source/style.rst index f26bcb85a..80fe8fd27 100644 --- a/docs/source/style.rst +++ b/docs/source/style.rst @@ -439,7 +439,7 @@ classes, or method definitions. Comment explaining example 1. - >>> np.add(1, 2) + >>> int(np.add(1, 2)) 3 Comment explaining a new example. diff --git a/examples/jnb.py b/examples/jnb.py index 5b8ead6ee..9135add1d 100644 --- a/examples/jnb.py +++ b/examples/jnb.py @@ -33,11 +33,12 @@ def py_file_to_string(src): if import_seen: # Once an import statement has been seen, break on encountering a line that # is neither an import statement nor a newline, nor a component of an import - # statement extended over multiple lines, nor an os.environ statement, nor - # components of a try/except construction (note that handling of these final - # two cases is probably not very robust). + # statement extended over multiple lines, nor an os.environ statement, nor a + # ray.init statement, nor components of a try/except construction (note that + # handling of these final two cases is probably not very robust). if not re.match( - r"(^import|^from|^\n$|^\W+[^\W]|^\)$|^os.environ|^try:$|^except)", line + r"(^import|^from|^\n$|^\W+[^\W]|^\)$|^os.environ|^ray.init|^try:$|^except)", + line, ): lines.append(line) break diff --git a/examples/notebooks_requirements.txt b/examples/notebooks_requirements.txt index bcb9e03c9..644a4db21 100644 --- a/examples/notebooks_requirements.txt +++ b/examples/notebooks_requirements.txt @@ -1,4 +1,6 @@ -r examples-requirements.txt +ipykernel +ipywidgets nbformat nbconvert nb_conda_kernels diff --git a/examples/scripts/README.rst b/examples/scripts/README.rst index 95d6aed8c..8c664c649 100644 --- a/examples/scripts/README.rst +++ b/examples/scripts/README.rst @@ -41,8 +41,6 @@ Computed Tomography CT Training and Reconstructions with UNet `ct_projector_comparison.py `_ X-ray Transform Comparison - `ct_multi_cs_tv_admm.py `_ - TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram) `ct_multi_tv_admm.py `_ TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors) @@ -176,6 +174,8 @@ Total Variation 3D TV-Regularized Sparse-View CT Reconstruction (Proximal ADMM Solver) `ct_astra_weighted_tv_admm.py `_ TV-Regularized Low-Dose CT Reconstruction + `ct_multi_tv_admm.py `_ + TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors) `ct_svmbir_tv_multi.py `_ TV-Regularized CT Reconstruction (Multiple Algorithms) `deconv_circ_tv_admm.py `_ @@ -275,6 +275,8 @@ ADMM 3D TV-Regularized Sparse-View CT Reconstruction (ADMM Solver) `ct_astra_weighted_tv_admm.py `_ TV-Regularized Low-Dose CT Reconstruction + `ct_multi_tv_admm.py `_ + TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors) `ct_svmbir_tv_multi.py `_ TV-Regularized CT Reconstruction (Multiple Algorithms) `ct_svmbir_ppp_bm3d_admm_cg.py `_ diff --git a/examples/scripts/ct_multi_cs_tv_admm.py b/examples/scripts/ct_multi_cs_tv_admm.py deleted file mode 100644 index cb8ddc6b6..000000000 --- a/examples/scripts/ct_multi_cs_tv_admm.py +++ /dev/null @@ -1,186 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# This file is part of the SCICO package. Details of the copyright -# and user license can be found in the 'LICENSE.txt' file distributed -# with the package. - -r""" -TV-Regularized Sparse-View CT Reconstruction (Multiple Projectors, Common Sinogram) -=================================================================================== - -This example demonstrates solution of a sparse-view CT reconstruction -problem with isotropic total variation (TV) regularization - - $$\mathrm{argmin}_{\mathbf{x}} \; (1/2) \| \mathbf{y} - A \mathbf{x} - \|_2^2 + \lambda \| C \mathbf{x} \|_{2,1} \;,$$ - -where $A$ is the X-ray transform (the CT forward projection operator), -$\mathbf{y}$ is the sinogram, $C$ is a 2D finite difference operator, and -$\mathbf{x}$ is the desired image. The solution is computed and compared -for all three 2D CT projectors available in scico, using a sinogram -computed with the svmbir projector. -""" - -import numpy as np - -from xdesign import Foam, discrete_phantom - -import scico.numpy as snp -from scico import functional, linop, loss, metric, plot -from scico.linop.xray import Parallel2dProjector, XRayTransform, astra, svmbir -from scico.optimize.admm import ADMM, LinearSubproblemSolver -from scico.util import device_info - -""" -Create a ground truth image. -""" -N = 512 # phantom size -np.random.seed(1234) -x_gt = snp.array(discrete_phantom(Foam(size_range=[0.075, 0.0025], gap=1e-3, porosity=1), size=N)) - -det_count = N -det_spacing = np.sqrt(2) - - -""" -Define CT geometry and construct array of (approximately) equivalent projectors. -""" -n_projection = 45 # number of projections -angles = np.linspace(0, np.pi, n_projection) # evenly spaced projection angles -projectors = { - "astra": astra.XRayTransform2D( - x_gt.shape, det_count, det_spacing, angles - np.pi / 2.0 - ), # astra - "svmbir": svmbir.XRayTransform( - x_gt.shape, 2 * np.pi - angles, det_count, delta_pixel=1.0, delta_channel=det_spacing - ), # svmbir - "scico": XRayTransform( - Parallel2dProjector((N, N), angles, det_count=det_count, dx=1 / det_spacing) - ), # scico -} - - -""" -Compute common sinogram using svmbir projector. -""" -A = projectors["astra"] -noise = np.random.normal(size=(n_projection, det_count)).astype(np.float32) -y = A @ x_gt + 2.0 * noise - - -""" -Construct initial solution for regularized problem. -""" -x0 = A.fbp(y) - - -""" -Solve the same problem using the different projectors. -""" -print(f"Solving on {device_info()}") -x_rec, hist = {}, {} -for p in projectors.keys(): - print(f"\nSolving with {p} projector") - - # Set up ADMM solver object. - λ = 2e1 # L1 norm regularization parameter - ρ = 1e3 # ADMM penalty parameter - maxiter = 100 # number of ADMM iterations - cg_tol = 1e-4 # CG relative tolerance - cg_maxiter = 50 # maximum CG iterations per ADMM iteration - - # The append=0 option makes the results of horizontal and vertical - # finite differences the same shape, which is required for the L21Norm, - # which is used so that g(Cx) corresponds to isotropic TV. - C = linop.FiniteDifference(input_shape=x_gt.shape, append=0) - g = λ * functional.L21Norm() - A = projectors[p] - f = loss.SquaredL2Loss(y=y, A=A) - - # Set up the solver. - solver = ADMM( - f=f, - g_list=[g], - C_list=[C], - rho_list=[ρ], - x0=x0, - maxiter=maxiter, - subproblem_solver=LinearSubproblemSolver(cg_kwargs={"tol": cg_tol, "maxiter": cg_maxiter}), - itstat_options={"display": True, "period": 5}, - ) - - # Run the solver. - solver.solve() - hist[p] = solver.itstat_object.history(transpose=True) - x_rec[p] = solver.x - - if p == "scico": - x_rec[p] = x_rec[p] * det_spacing # to match ASTRA's scaling - - -""" -Compare reconstruction results. -""" -print("Reconstruction SNR:") -for p in projectors.keys(): - print(f" {(p + ':'):7s} {metric.snr(x_gt, x_rec[p]):5.2f} dB") - - -""" -Display sinogram. -""" -fig, ax = plot.subplots(nrows=1, ncols=1, figsize=(15, 3)) -plot.imview(y, title="sinogram", fig=fig, ax=ax) -fig.show() - - -""" -Plot convergence statistics. -""" -fig, ax = plot.subplots(nrows=1, ncols=3, figsize=(12, 5)) -plot.plot( - np.vstack([hist[p].Objective for p in projectors.keys()]).T, - title="Objective function", - xlbl="Iteration", - ylbl="Functional value", - lgnd=projectors.keys(), - fig=fig, - ax=ax[0], -) -plot.plot( - np.vstack([hist[p].Prml_Rsdl for p in projectors.keys()]).T, - ptyp="semilogy", - title="Primal Residual", - xlbl="Iteration", - fig=fig, - ax=ax[1], -) -plot.plot( - np.vstack([hist[p].Dual_Rsdl for p in projectors.keys()]).T, - ptyp="semilogy", - title="Dual Residual", - xlbl="Iteration", - fig=fig, - ax=ax[2], -) -fig.show() - - -""" -Show the recovered images. -""" -fig, ax = plot.subplots(nrows=1, ncols=4, figsize=(15, 5)) -plot.imview(x_gt, title="Ground truth", fig=fig, ax=ax[0]) -for n, p in enumerate(projectors.keys()): - plot.imview( - x_rec[p], - title="%s SNR: %.2f (dB)" % (p, metric.snr(x_gt, x_rec[p])), - fig=fig, - ax=ax[n + 1], - ) -for ax in ax: - ax.get_images()[0].set_clim(-0.1, 1.1) -fig.show() - - -input("\nWaiting for input to close figures and exit") diff --git a/examples/scripts/ct_multi_tv_admm.py b/examples/scripts/ct_multi_tv_admm.py index 084ad99ac..87fc86865 100644 --- a/examples/scripts/ct_multi_tv_admm.py +++ b/examples/scripts/ct_multi_tv_admm.py @@ -17,7 +17,8 @@ where $A$ is the X-ray transform (the CT forward projection operator), $\mathbf{y}$ is the sinogram, $C$ is a 2D finite difference operator, and $\mathbf{x}$ is the desired image. The solution is computed and compared -for all three 2D CT projectors available in scico. +for all three 2D CT projectors available in scico, using a sinogram +computed with the astra projector. """ import numpy as np @@ -37,6 +38,9 @@ np.random.seed(1234) x_gt = snp.array(discrete_phantom(Foam(size_range=[0.075, 0.0025], gap=1e-3, porosity=1), size=N)) +det_count = N +det_spacing = np.sqrt(2) + """ Define CT geometry and construct array of (approximately) equivalent projectors. @@ -44,37 +48,54 @@ n_projection = 45 # number of projections angles = np.linspace(0, np.pi, n_projection) # evenly spaced projection angles projectors = { - "astra": astra.XRayTransform2D(x_gt.shape, N, 1.0, angles - np.pi / 2.0), # astra - "svmbir": svmbir.XRayTransform(x_gt.shape, 2 * np.pi - angles, N), # svmbir - "scico": XRayTransform(Parallel2dProjector((N, N), angles, det_count=N)), # scico + "astra": astra.XRayTransform2D( + x_gt.shape, det_count, det_spacing, angles - np.pi / 2.0 + ), # astra + "svmbir": svmbir.XRayTransform( + x_gt.shape, 2 * np.pi - angles, det_count, delta_pixel=1.0, delta_channel=det_spacing + ), # svmbir + "scico": XRayTransform( + Parallel2dProjector((N, N), angles, det_count=det_count, dx=1 / det_spacing) + ), # scico } +""" +Compute common sinogram using astra projector. +""" +A = projectors["astra"] +noise = np.random.normal(size=(n_projection, det_count)).astype(np.float32) +y = A @ x_gt + 2.0 * noise + + +""" +Construct initial solution for regularized problem. +""" +x0 = A.fbp(y) + + """ Solve the same problem using the different projectors. """ print(f"Solving on {device_info()}") -y, x_rec, hist = {}, {}, {} -noise = np.random.normal(size=(n_projection, N)).astype(np.float32) -for p in ("astra", "svmbir", "scico"): +x_rec, hist = {}, {} +for p in projectors.keys(): print(f"\nSolving with {p} projector") - A = projectors[p] - y[p] = A @ x_gt + 2.0 * noise # sinogram # Set up ADMM solver object. - λ = 2e0 # L1 norm regularization parameter - ρ = 5e0 # ADMM penalty parameter - maxiter = 25 # number of ADMM iterations + λ = 2e1 # L1 norm regularization parameter + ρ = 1e3 # ADMM penalty parameter + maxiter = 100 # number of ADMM iterations cg_tol = 1e-4 # CG relative tolerance - cg_maxiter = 25 # maximum CG iterations per ADMM iteration + cg_maxiter = 50 # maximum CG iterations per ADMM iteration # The append=0 option makes the results of horizontal and vertical # finite differences the same shape, which is required for the L21Norm, # which is used so that g(Cx) corresponds to isotropic TV. C = linop.FiniteDifference(input_shape=x_gt.shape, append=0) g = λ * functional.L21Norm() - f = loss.SquaredL2Loss(y=y[p], A=A) - x0 = snp.clip(A.T(y[p]), 0, 1.0) + A = projectors[p] + f = loss.SquaredL2Loss(y=y, A=A) # Set up the solver. solver = ADMM( @@ -91,15 +112,25 @@ # Run the solver. solver.solve() hist[p] = solver.itstat_object.history(transpose=True) - x_rec[p] = snp.clip(solver.x, 0, 1.0) + x_rec[p] = solver.x + + if p == "scico": + x_rec[p] = x_rec[p] * det_spacing # to match ASTRA's scaling + + +""" +Compare reconstruction results. +""" +print("Reconstruction SNR:") +for p in projectors.keys(): + print(f" {(p + ':'):7s} {metric.snr(x_gt, x_rec[p]):5.2f} dB") """ -Compare sinograms. +Display sinogram. """ -fig, ax = plot.subplots(nrows=3, ncols=1, figsize=(15, 10)) -for idx, name in enumerate(projectors.keys()): - plot.imview(y[name], title=f"{name} sinogram", cbar=None, fig=fig, ax=ax[idx]) +fig, ax = plot.subplots(nrows=1, ncols=1, figsize=(15, 3)) +plot.imview(y, title="sinogram", fig=fig, ax=ax) fig.show() @@ -147,6 +178,8 @@ fig=fig, ax=ax[n + 1], ) +for ax in ax: + ax.get_images()[0].set_clim(-0.1, 1.1) fig.show() diff --git a/examples/scripts/deconv_circ_tv_admm.py b/examples/scripts/deconv_circ_tv_admm.py index b2ba83202..90e898dc7 100644 --- a/examples/scripts/deconv_circ_tv_admm.py +++ b/examples/scripts/deconv_circ_tv_admm.py @@ -54,7 +54,7 @@ """ Set up an ADMM solver object. """ -λ = 2e-2 # L21 norm regularization parameter +λ = 2e-2 # ℓ2,1 norm regularization parameter ρ = 5e-1 # ADMM penalty parameter maxiter = 50 # number of ADMM iterations diff --git a/examples/scripts/deconv_microscopy_allchn_tv_admm.py b/examples/scripts/deconv_microscopy_allchn_tv_admm.py index 13ec27b87..ccd8c9471 100644 --- a/examples/scripts/deconv_microscopy_allchn_tv_admm.py +++ b/examples/scripts/deconv_microscopy_allchn_tv_admm.py @@ -28,19 +28,24 @@ non-negativity constraint, and $\mathbf{x}$ is the desired image. """ +# isort: off import numpy as np +import logging import ray + +ray.init(logging_level=logging.ERROR) # need to call init before jax import: ray-project/ray#44087 + import scico.numpy as snp from scico import functional, linop, loss, plot from scico.examples import downsample_volume, epfl_deconv_data, tile_volume_slices from scico.optimize.admm import ADMM, CircularConvolveSolver """ -Get and preprocess data. We downsample the data for the for purposes of -the example. Reducing the downsampling rate will make the example slower -and more memory-intensive. To run this example on a GPU it may be -necessary to set environment variables +Get and preprocess data. The data is downsampled to limit the memory +requirements and run time of the example. Reducing the downsampling rate +will make the example slower and more memory-intensive. To run this +example on a GPU it may be necessary to set environment variables `XLA_PYTHON_CLIENT_ALLOCATOR=platform` and `XLA_PYTHON_CLIENT_PREALLOCATE=false`. If your GPU does not have enough memory, you can try setting the environment variable @@ -81,11 +86,9 @@ """ -Initialize ray, determine available computing resources, and put large arrays -in object store. +Determine available computing resources, and put large arrays in ray +object store. """ -ray.init() - ngpu = 0 ar = ray.available_resources() ncpu = max(int(ar["CPU"]) // 3, 1) diff --git a/examples/scripts/deconv_microscopy_tv_admm.py b/examples/scripts/deconv_microscopy_tv_admm.py index eee4d8f42..5c9180b60 100644 --- a/examples/scripts/deconv_microscopy_tv_admm.py +++ b/examples/scripts/deconv_microscopy_tv_admm.py @@ -34,10 +34,10 @@ from scico.optimize.admm import ADMM, CircularConvolveSolver """ -Get and preprocess data. We downsample the data for the for purposes of -the example. Reducing the downsampling rate will make the example slower -and more memory-intensive. To run this example on a GPU it may be -necessary to set environment variables +Get and preprocess data. The data is downsampled to limit the memory +requirements and run time of the example. Reducing the downsampling rate +will make the example slower and more memory-intensive. To run this +example on a GPU it may be necessary to set environment variables `XLA_PYTHON_CLIENT_ALLOCATOR=platform` and `XLA_PYTHON_CLIENT_PREALLOCATE=false`. If your GPU does not have enough memory, you can try setting the environment variable diff --git a/examples/scripts/deconv_tv_admm.py b/examples/scripts/deconv_tv_admm.py index 874b70e3c..8acb54d36 100644 --- a/examples/scripts/deconv_tv_admm.py +++ b/examples/scripts/deconv_tv_admm.py @@ -77,7 +77,7 @@ f = loss.SquaredL2Loss(y=y, A=C) # Penalty parameters must be accounted for in the gi functions, not as # additional inputs. -λ = 2.1e-2 # L21 norm regularization parameter +λ = 2.1e-2 # ℓ2,1 norm regularization parameter g = λ * functional.L21Norm() # The append=0 option makes the results of horizontal and vertical # finite differences the same shape, which is required for the L21Norm, diff --git a/examples/scripts/index.rst b/examples/scripts/index.rst index bbc697a95..8100f06e6 100644 --- a/examples/scripts/index.rst +++ b/examples/scripts/index.rst @@ -25,7 +25,6 @@ Computed Tomography - ct_astra_odp_train_foam2.py - ct_astra_unet_train_foam2.py - ct_projector_comparison.py - - ct_multi_cs_tv_admm.py - ct_multi_tv_admm.py Deconvolution @@ -106,6 +105,7 @@ Total Variation - ct_astra_3d_tv_admm.py - ct_astra_3d_tv_padmm.py - ct_astra_weighted_tv_admm.py + - ct_multi_tv_admm.py - ct_svmbir_tv_multi.py - deconv_circ_tv_admm.py - deconv_tv_admm.py @@ -165,6 +165,7 @@ ADMM - ct_tv_admm.py - ct_astra_3d_tv_admm.py - ct_astra_weighted_tv_admm.py + - ct_multi_tv_admm.py - ct_svmbir_tv_multi.py - ct_svmbir_ppp_bm3d_admm_cg.py - ct_svmbir_ppp_bm3d_admm_prox.py diff --git a/examples/scripts/sparsecode_conv_admm.py b/examples/scripts/sparsecode_conv_admm.py index 611e7238b..0c46356a3 100644 --- a/examples/scripts/sparsecode_conv_admm.py +++ b/examples/scripts/sparsecode_conv_admm.py @@ -75,7 +75,7 @@ """ Set functional and solver parameters. """ -λ = 1e0 # l1-l2 norm regularization parameter +λ = 1e0 # ℓ1-ℓ2 norm regularization parameter ρ = 2e0 # ADMM penalty parameter maxiter = 200 # number of ADMM iterations diff --git a/examples/scripts/sparsecode_conv_md_admm.py b/examples/scripts/sparsecode_conv_md_admm.py index 9f70e9ce3..77d24f19b 100644 --- a/examples/scripts/sparsecode_conv_md_admm.py +++ b/examples/scripts/sparsecode_conv_md_admm.py @@ -96,7 +96,7 @@ """ Set functional and solver parameters. """ -λ = 1e0 # l1-l2 norm regularization parameter +λ = 1e0 # ℓ1-ℓ2 norm regularization parameter ρ0 = 1e0 # ADMM penalty parameters ρ1 = 3e0 maxiter = 200 # number of ADMM iterations diff --git a/examples/scripts/video_rpca_admm.py b/examples/scripts/video_rpca_admm.py index e6c71cb8e..ee8fa5c14 100644 --- a/examples/scripts/video_rpca_admm.py +++ b/examples/scripts/video_rpca_admm.py @@ -62,7 +62,7 @@ Set up an ADMM solver object. """ λ0 = 1e1 # nuclear norm regularization parameter -λ1 = 3e1 # l1 norm regularization parameter +λ1 = 3e1 # ℓ1 norm regularization parameter ρ0 = 2e1 # ADMM penalty parameter ρ1 = 2e1 # ADMM penalty parameter maxiter = 50 # number of ADMM iterations diff --git a/scico/optimize/_pgm.py b/scico/optimize/_pgm.py index 8375c0b2d..6188f7480 100644 --- a/scico/optimize/_pgm.py +++ b/scico/optimize/_pgm.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright (C) 2020-2023 by SCICO Developers +# Copyright (C) 2020-2024 by SCICO Developers # All rights reserved. BSD 3-clause License. # This file is part of the SCICO package. Details of the copyright and # user license can be found in the 'LICENSE' file distributed with the @@ -30,14 +30,35 @@ class PGM(Optimizer): - r"""Proximal Gradient Method (PGM) base class. + r"""Proximal gradient method (PGM) algorithm. - Minimize a function of the form :math:`f(\mb{x}) + g(\mb{x})`, where - :math:`f` and the :math:`g` are instances of :class:`.Functional`. + Minimize a functional of the form :math:`f(\mb{x}) + g(\mb{x})`, + where :math:`f` and the :math:`g` are instances of + :class:`.Functional`. Functional :math:`f` should be differentiable + and have a Lipschitz continuous derivative, and functional :math:`g` + should have a proximal operator defined. + + The step size :math:`\alpha` of the algorithm is defined in terms of + its reciprocal :math:`L`, i.e. :math:`\alpha = 1 / L`. The initial + value for this parameter, `L0`, is required to satisfy + + .. math:: + L_0 \geq K(\nabla f) \;, - Uses helper :class:`StepSize` to provide an estimate of the Lipschitz - constant :math:`L` of :math:`f`. The step size :math:`\alpha` is the - reciprocal of :math:`L`, i.e.: :math:`\alpha = 1 / L`. + where :math:`K(\nabla f)` denotes the Lipschitz constant of the + gradient of :math:`f`. When `f` is an instance of + :class:`.SquaredL2Loss` with a :class:`.LinearOperator` `A`, + + .. math:: + K(\nabla f) = \lambda_{ \mathrm{max} }( A^H A ) = \| A \|_2^2 \;, + + where :math:`\lambda_{\mathrm{max}}(B)` denotes the largest + eigenvalue of :math:`B`. + + The evolution of the step size is controlled by auxiliary class + :class:`.PGMStepSize` and derived classes. The default + :class:`.PGMStepSize` simply sets :math:`L = L_0`, while the derived + classes implement a variety of adaptive strategies. """ def __init__( @@ -52,12 +73,14 @@ def __init__( r""" Args: - f: Loss or Functional object with `grad` defined. - g: Instance of Functional with defined prox method. - L0: Initial estimate of Lipschitz constant of f. + f: Instance of :class:`.Loss` or :class:`.Functional` with + defined `grad` method. + g: Instance of :class:`.Functional` with defined prox method. + L0: Initial estimate of Lipschitz constant of gradient of `f`. x0: Starting point for :math:`\mb{x}`. - step_size: helper :class:`StepSize` to estimate the Lipschitz - constant of f. + step_size: Instance of an auxiliary class of type + :class:`.PGMStepSize` determining the evolution of the + algorithm step size. **kwargs: Additional optional parameters handled by initializer of base class :class:`.Optimizer`. """ @@ -75,7 +98,7 @@ def __init__( step_size = PGMStepSize() self.step_size: PGMStepSize = step_size self.step_size.internal_init(self) - self.L: float = L0 # reciprocal of step size (estimate of Lipschitz constant of f) + self.L: float = L0 # reciprocal of step size (estimate of Lipschitz constant of ∇f) self.fixed_point_residual = snp.inf def x_step(v: Union[Array, BlockArray], L: float) -> Union[Array, BlockArray]: @@ -150,16 +173,14 @@ def step(self): class AcceleratedPGM(PGM): - r"""Accelerated Proximal Gradient Method (AcceleratedPGM) base class. - - Minimize a function of the form :math:`f(\mb{x}) + g(\mb{x})`. + r"""Accelerated proximal gradient method (APGM) algorithm. Minimize a function of the form :math:`f(\mb{x}) + g(\mb{x})`, where :math:`f` and the :math:`g` are instances of :class:`.Functional`. The accelerated form of PGM is also known as FISTA :cite:`beck-2009-fast`. - For documentation on inherited attributes, see :class:`.PGM`. + See :class:`.PGM` for more detailed documentation. """ def __init__( @@ -173,12 +194,14 @@ def __init__( ): r""" Args: - f: Loss or Functional object with `grad` defined. - g: Instance of Functional with defined prox method. - L0: Initial estimate of Lipschitz constant of f. + f: Instance of :class:`.Loss` or :class:`.Functional` with + defined `grad` method. + g: Instance of :class:`.Functional` with defined prox method. + L0: Initial estimate of Lipschitz constant of gradient of `f`. x0: Starting point for :math:`\mb{x}`. - step_size: helper :class:`StepSize` to estimate the Lipschitz - constant of f. + step_size: Instance of an auxiliary class of type + :class:`.PGMStepSize` determining the evolution of the + algorithm step size. **kwargs: Additional optional parameters handled by initializer of base class :class:`.Optimizer`. """ diff --git a/scico/plot.py b/scico/plot.py index 1cdd1adc7..6c0375be1 100644 --- a/scico/plot.py +++ b/scico/plot.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright (C) 2020-2022 by SCICO Developers +# Copyright (C) 2020-2024 by SCICO Developers # All rights reserved. BSD 3-clause License. # This file is part of the SCICO package. Details of the copyright and # user license can be found in the 'LICENSE' file distributed with the @@ -360,7 +360,7 @@ def surf( # https://stackoverflow.com/a/35221116 if ax.name != "3d": ax.remove() - ax = fig.add_subplot(*ax.get_geometry(), projection="3d") + ax = fig.add_subplot(ax.get_subplotspec(), projection="3d") if elev is not None or azim is not None: ax.view_init(elev=elev, azim=azim)