Convergence for nonconvex ADMM, with applications to CT imaging.

The alternating direction method of multipliers (ADMM) algorithm is a powerful and flexible tool for complex optimization problems of the form m i n { f ( x ) + g ( y ) : A x + B y = c } . ADMM exhibits robust empirical performance across a range of challenging settings including nonsmoothness and nonconvexity of the objective functions f and g , and provides a simple and natural approach to the inverse problem of image reconstruction for computed tomography (CT) imaging. From the theoretical point of view, existing results for convergence in the nonconvex setting generally assume smoothness in at least one of the component functions in the objective. In this work, our new theoretical results provide convergence guarantees under a restricted strong convexity assumption without requiring smoothness or differentiability, while still allowing differentiable terms to be treated approximately if needed. We validate these theoretical results empirically, with a simulated example where both f and g are nondifferentiable-and thus outside the scope of existing theory-as well as a simulated CT image reconstruction problem.

Simultaneous Activity and Attenuation Estimation in TOF-PET With TV-Constrained Nonconvex Optimization.

An alternating direction method of multipliers (ADMM) framework is developed for nonsmooth biconvex optimization for inverse problems in imaging. In particular, the simultaneous estimation of activity and attenuation (SAA) problem in time-of-flight positron emission tomography (TOF-PET) has such a structure when maximum likelihood estimation (MLE) is employed. The ADMM framework is applied to MLE for SAA in TOF-PET, resulting in the ADMM-SAA algorithm. This algorithm is extended by imposing total variation (TV) constraints on both the activity and attenuation map, resulting in the ADMM-TVSAA algorithm. The performance of this algorithm is illustrated using the penalized maximum likelihood activity and attenuation estimation (P-MLAA) algorithm as a reference.

Simultaneous activity and attenuation estimation in TOF-PET with TV-constrained nonconvex optimization.

An alternating direction method of multipliers (ADMM) framework is developed for nonsmooth biconvex optimization for inverse problems in imaging. In particular, the simultaneous estimation of activity and attenuation (SAA) problem in time-of-flight positron emission tomography (TOF-PET) has such a structure when maximum likelihood estimation (MLE) is employed. The ADMM framework is applied to MLE for SAA in TOF-PET, resulting in the ADMM-SAA algorithm. This algorithm is extended by imposing total variation (TV) constraints on both the activity and attenuation map, resulting in the ADMM-TVSAA algorithm. The performance of this algorithm is illustrated using the penalized maximum likelihood activity and attenuation estimation (P-MLAA) algorithm as a reference. Additional results on step-size tuning and on the use of unconstrained ADMM-SAA are presented in the previous arXiv submission: arXiv:2303.17042v1.

Black-box tests for algorithmic stability.

Algorithmic stability is a concept from learning theory that expresses the degree to which changes to the input data (e.g. removal of a single data point) may affect the outputs of a regression algorithm. Knowing an algorithm's stability properties is often useful for many downstream applications-for example, stability is known to lead to desirable generalization properties and predictive inference guarantees. However, many modern algorithms currently used in practice are too complex for a theoretical analysis of their stability properties, and thus we can only attempt to establish these properties through an empirical exploration of the algorithm's behaviour on various datasets. In this work, we lay out a formal statistical framework for this kind of black-box testing without any assumptions on the algorithm or the data distribution, and establish fundamental bounds on the ability of any black-box test to identify algorithmic stability.

Constrained one-step material decomposition reconstruction of head CT data from a silicon photon-counting prototype.

BACKGROUND: Spectral CT material decomposition provides quantitative information but is challenged by the instability of the inversion into basis materials. We have previously proposed the constrained One-Step Spectral CT Image Reconstruction (cOSSCIR) algorithm to stabilize the material decomposition inversion by directly estimating basis material images from spectral CT data. cOSSCIR was previously investigated on phantom data. PURPOSE: This study investigates the performance of cOSSCIR using head CT datasets acquired on a clinical photon-counting CT (PCCT) prototype. This is the first investigation of cOSSCIR for large-scale, anatomically complex, clinical PCCT data. The cOSSCIR decomposition is preceded by a spectrum estimation and nonlinear counts correction calibration step to address nonideal detector effects. METHODS: Head CT data were acquired on an early prototype clinical PCCT system using an edge-on silicon detector with eight energy bins. Calibration data of a step wedge phantom were also acquired and used to train a spectral model to account for the source spectrum and detector spectral response, and also to train a nonlinear counts correction model to account for pulse pileup effects. The cOSSCIR algorithm optimized the bone and adipose basis images directly from the photon counts data, while placing a grouped total variation (TV) constraint on the basis images. For comparison, basis images were also reconstructed by a two-step projection-domain approach of Maximum Likelihood Estimation (MLE) for decomposing basis sinograms, followed by filtered backprojection (MLE + FBP) or a TV minimization algorithm (MLE + TVmin ) to reconstruct basis images. We hypothesize that the cOSSCIR approach will provide a more stable inversion into basis images compared to two-step approaches. To investigate this hypothesis, the noise standard deviation in bone and soft-tissue regions of interest (ROIs) in the reconstructed images were compared between cOSSCIR and the two-step methods for a range of regularization constraint settings. RESULTS: cOSSCIR reduced the noise standard deviation in the basis images by a factor of two to six compared to that of MLE + TVmin , when both algorithms were constrained to produce images with the same TV. The cOSSCIR images demonstrated qualitatively improved spatial resolution and depiction of fine anatomical detail. The MLE + TVmin algorithm resulted in lower noise standard deviation than cOSSCIR for the virtual monoenergetic images (VMIs) at higher energy levels and constraint settings, while the cOSSCIR VMIs resulted in lower noise standard deviation at lower energy levels and overall higher qualitative spatial resolution. There were no statistically significant differences in the mean values within the bone region of images reconstructed by the studied algorithms. There were statistically significant differences in the mean values within the soft-tissue region of the reconstructed images, with cOSSCIR producing mean values closer to the expected values. CONCLUSIONS: The cOSSCIR algorithm, combined with our previously proposed spectral model estimation and nonlinear counts correction method, successfully estimated bone and adipose basis images from high resolution, large-scale patient data from a clinical PCCT prototype. The cOSSCIR basis images were able to depict fine anatomical details with a factor of two to six reduction in noise standard deviation compared to that of the MLE + TVmin two-step approach.

Addressing CT metal artifacts using photon-counting detectors and one-step spectral CT image reconstruction.

PURPOSE: The constrained one-step spectral CT image reconstruction (cOSSCIR) algorithm with a nonconvex alternating direction method of multipliers optimizer is proposed for addressing computed tomography (CT) metal artifacts caused by beam hardening, noise, and photon starvation. The quantitative performance of cOSSCIR is investigated through a series of photon-counting CT simulations. METHODS: cOSSCIR directly estimates basis material maps from photon-counting data using a physics-based forward model that accounts for beam hardening. The cOSSCIR optimization framework places constraints on the basis maps, which we hypothesize will stabilize the decomposition and reduce streaks caused by noise and photon starvation. Another advantage of cOSSCIR is that the spectral data need not be registered, so that a ray can be used even if some energy window measurements are unavailable. Photon-counting CT acquisitions of a virtual pelvic phantom with low-contrast soft tissue texture and bilateral hip prostheses were simulated. Bone and water basis maps were estimated using the cOSSCIR algorithm and combined to form a virtual monoenergetic image for the evaluation of metal artifacts. The cOSSCIR images were compared to a "two-step" decomposition approach that first estimated basis sinograms using a maximum likelihood algorithm and then reconstructed basis maps using an iterative total variation constrained least-squares optimization (MLE+TV min $_{\text{min}}$ ). Images were also compared to a nonspectral TV min $_{\text{min}}$ reconstruction of the total number of counts detected for each ray with and without normalized metal artifact reduction (NMAR) applied. The simulated metal density was increased to investigate the effects of increasing photon starvation. The quantitative error and standard deviation in regions of the phantom were compared across the investigated algorithms. The ability of cOSSCIR to reproduce the soft-tissue texture, while reducing metal artifacts, was quantitatively evaluated. RESULTS: Noiseless simulations demonstrated the convergence of the cOSSCIR and MLE+TV min $_{\text{min}}$ algorithms to the correct basis maps in the presence of beam-hardening effects. When noise was simulated, cOSSCIR demonstrated a quantitative error of -1 HU, compared to 2 HU error for the MLE+TV min $_{\text{min}}$ algorithm and -154 HU error for the nonspectral TV min $_{\text{min}}$ +NMAR algorithm. For the cOSSCIR algorithm, the standard deviation in the central iodine region of interest was 20 HU, compared to 299 HU for the MLE+TV min $_{\text{min}}$ algorithm, 41 HU for the MLE+TV min $_{\text{min}}$ +Mask algorithm that excluded rays through metal, and 55 HU for the nonspectral TV min $_{\text{min}}$ +NMAR algorithm. Increasing levels of photon starvation did not impact the bias or standard deviation of the cOSSCIR images. cOSSCIR was able to reproduce the soft-tissue texture when an appropriate regularization constraint value was selected. CONCLUSIONS: By directly inverting photon-counting CT data into basis maps using an accurate physics-based forward model and a constrained optimization algorithm, cOSSCIR avoids metal artifacts due to beam hardening, noise, and photon starvation. The cOSSCIR algorithm demonstrated improved stability and accuracy compared to a two-step method of decomposition followed by reconstruction.

Fast and flexible estimation of effective migration surfaces.

Spatial population genetic data often exhibits 'isolation-by-distance,' where genetic similarity tends to decrease as individuals become more geographically distant. The rate at which genetic similarity decays with distance is often spatially heterogeneous due to variable population processes like genetic drift, gene flow, and natural selection. Petkova et al., 2016 developed a statistical method called Estimating Effective Migration Surfaces (EEMS) for visualizing spatially heterogeneous isolation-by-distance on a geographic map. While EEMS is a powerful tool for depicting spatial population structure, it can suffer from slow runtimes. Here, we develop a related method called Fast Estimation of Effective Migration Surfaces (FEEMS). FEEMS uses a Gaussian Markov Random Field model in a penalized likelihood framework that allows for efficient optimization and output of effective migration surfaces. Further, the efficient optimization facilitates the inference of migration parameters per edge in the graph, rather than per node (as in EEMS). With simulations, we show conditions under which FEEMS can accurately recover effective migration surfaces with complex gene-flow histories, including those with anisotropy. We apply FEEMS to population genetic data from North American gray wolves and show it performs favorably in comparison to EEMS, with solutions obtained orders of magnitude faster. Overall, FEEMS expands the ability of users to quickly visualize and interpret spatial structure in their data.

The bias of isotonic regression.

We study the bias of the isotonic regression estimator. While there is extensive work characterizing the mean squared error of the isotonic regression estimator, relatively little is known about the bias. In this paper, we provide a sharp characterization, proving that the bias scales as O(n -ß/3) up to log factors, where 1 ≤ ß ≤ 2 is the exponent corresponding to Hölder smoothness of the underlying mean. Importantly, this result only requires a strictly monotone mean and that the noise distribution has subexponential tails, without relying on symmetric noise or other restrictive assumptions.

Estimating the spectrum in computed tomography via Kullback-Leibler divergence constrained optimization.

PURPOSE: We study the problem of spectrum estimation from transmission data of a known phantom. The goal is to reconstruct an x-ray spectrum that can accurately model the x-ray transmission curves and reflects a realistic shape of the typical energy spectra of the CT system. METHODS: Spectrum estimation is posed as an optimization problem with x-ray spectrum as unknown variables, and a Kullback-Leibler (KL)-divergence constraint is employed to incorporate prior knowledge of the spectrum and enhance numerical stability of the estimation process. The formulated constrained optimization problem is convex and can be solved efficiently by use of the exponentiated-gradient (EG) algorithm. We demonstrate the effectiveness of the proposed approach on the simulated and experimental data. The comparison to the expectation-maximization (EM) method is also discussed. RESULTS: In simulations, the proposed algorithm is seen to yield x-ray spectra that closely match the ground truth and represent the attenuation process of x-ray photons in materials, both included and not included in the estimation process. In experiments, the calculated transmission curve is in good agreement with the measured transmission curve, and the estimated spectra exhibits physically realistic looking shapes. The results further show the comparable performance between the proposed optimization-based approach and EM. CONCLUSIONS: Our formulation of a constrained optimization provides an interpretable and flexible framework for spectrum estimation. Moreover, a KL-divergence constraint can include a prior spectrum and appears to capture important features of x-ray spectrum, allowing accurate and robust estimation of x-ray spectrum in CT imaging.

EigenPrism: inference for high dimensional signal-to-noise ratios.

Consider the following three important problems in statistical inference, namely, constructing confidence intervals for (1) the error of a high-dimensional (p > n) regression estimator, (2) the linear regression noise level, and (3) the genetic signal-to-noise ratio of a continuous-valued trait (related to the heritability). All three problems turn out to be closely related to the little-studied problem of performing inference on the [Formula: see text]-norm of the signal in high-dimensional linear regression. We derive a novel procedure for this, which is asymptotically correct when the covariates are multivariate Gaussian and produces valid confidence intervals in finite samples as well. The procedure, called EigenPrism, is computationally fast and makes no assumptions on coefficient sparsity or knowledge of the noise level. We investigate the width of the EigenPrism confidence intervals, including a comparison with a Bayesian setting in which our interval is just 5% wider than the Bayes credible interval. We are then able to unify the three aforementioned problems by showing that the EigenPrism procedure with only minor modifications is able to make important contributions to all three. We also investigate the robustness of coverage and find that the method applies in practice and in finite samples much more widely than just the case of multivariate Gaussian covariates. Finally, we apply EigenPrism to a genetic dataset to estimate the genetic signal-to-noise ratio for a number of continuous phenotypes.

A Spectral CT Method to Directly Estimate Basis Material Maps From Experimental Photon-Counting Data.

The proposed spectral CT method solves the constrained one-step spectral CT reconstruction (cOSSCIR) optimization problem to estimate basis material maps while modeling the nonlinear X-ray detection process and enforcing convex constraints on the basis map images. In order to apply the optimization-based reconstruction approach to experimental data, the presented method empirically estimates the effective energy-window spectra using a calibration procedure. The amplitudes of the estimated spectra were further optimized as part of the reconstruction process to reduce ring artifacts. A validation approach was developed to select constraint parameters. The proposed spectral CT method was evaluated through simulations and experiments with a photon-counting detector. Basis material map images were successfully reconstructed using the presented empirical spectral modeling and cOSSCIR optimization approach. In simulations, the cOSSCIR approach accurately reconstructed the basis map images (<1% error). In experiments, the proposed method estimated the low-density polyethylene region of the basis maps with 0.5% error in the PMMA image and 4% error in the aluminum image. For the Teflon region, the experimental results demonstrated 8% and 31% error in the PMMA and aluminum basis material maps, respectively, compared with -24% and 126% error without estimation of the effective energy window spectra, with residual errors likely due to insufficient modeling of detector effects. The cOSSCIR algorithm estimated the material decomposition angle to within 1.3 degree error, where, for reference, the difference in angle between PMMA and muscle tissue is 2.1 degrees. The joint estimation of spectral-response scaling coefficients and basis material maps was found to reduce ring artifacts in both a phantom and tissue specimen. The presented validation procedure demonstrated feasibility for the automated determination of algorithm constraint parameters.

MOCCA: Mirrored Convex/Concave Optimization for Nonconvex Composite Functions.

Many optimization problems arising in high-dimensional statistics decompose naturally into a sum of several terms, where the individual terms are relatively simple but the composite objective function can only be optimized with iterative algorithms. In this paper, we are interested in optimization problems of the form F(Kx) + G(x), where K is a fixed linear transformation, while F and G are functions that may be nonconvex and/or nondifferentiable. In particular, if either of the terms are nonconvex, existing alternating minimization techniques may fail to converge; other types of existing approaches may instead be unable to handle nondifferentiability. We propose the MOCCA (mirrored convex/concave) algorithm, a primal/dual optimization approach that takes a local convex approximation to each term at every iteration. Inspired by optimization problems arising in computed tomography (CT) imaging, this algorithm can handle a range of nonconvex composite optimization problems, and offers theoretical guarantees for convergence when the overall problem is approximately convex (that is, any concavity in one term is balanced out by convexity in the other term). Empirical results show fast convergence for several structured signal recovery problems.