Misha E. Kilmer

Meghan O'Connell, Misha E. Kilmer, Eric de Sturler et al.

In nonlinear imaging problems whose forward model is described by a partial differential equation (PDE), the main computational bottleneck in solving the inverse problem is the need to solve many large-scale discretized PDEs at each step of the optimization process. In the context of absorption imaging in diffuse optical tomography, one approach to addressing this bottleneck proposed recently (de Sturler, et al, 2015) reformulates the viewing of the forward problem as a differential algebraic system, and then employs model order reduction (MOR). However, the construction of the reduced model requires the solution of several full order problems (i.e. the full discretized PDE for multiple right-hand sides) to generate a candidate global basis. This step is then followed by a rank-revealing factorization of the matrix containing the candidate basis in order to compress the basis to a size suitable for constructing the reduced transfer function. The present paper addresses the costs associated with the global basis approximation in two ways. First, we use the structure of the matrix to rewrite the full order transfer function, and corresponding derivatives, such that the full order systems to be solved are symmetric (positive definite in the zero frequency case). Then we apply MOR to the new formulation of the problem. Second, we give an approach to computing the global basis approximation dynamically as the full order systems are solved. In this phase, only the incrementally new, relevant information is added to the existing global basis, and redundant information is not computed. This new approach is achieved by an inner-outer Krylov recycling approach which has potential use in other applications as well. We show the value of the new approach to approximate global basis computation on two DOT absorption image reconstruction problems.

Selin Aslan, Eric de Sturler, Misha E. Kilmer

In partial differential equations-based (PDE-based) inverse problems with many measurements, many large-scale discretized PDEs must be solved for each evaluation of the misfit or objective function. In the nonlinear case, evaluating the Jacobian requires solving an additional set of systems. This leads to a tremendous computational cost, and this is by far the dominant cost for these problems. Several authors have proposed randomization and stochastic programming techniques to drastically reduce the number of system solves by estimating the objective function using only a few appropriately chosen random linear combinations of the sources. While some have reported good solution quality at a greatly reduced cost, for our problem of interest, diffuse optical tomography, the approach often does not lead to sufficiently accurate solutions. We propose two improvements. First, to efficiently exploit Newton-type methods, we modify the stochastic estimates to include random linear combinations of detectors, drastically reducing the number of adjoint solves. Second, after solving to a modest tolerance, we compute a few simultaneous sources and detectors that maximize the Frobenius norm of the sampled Jacobian to improve the rate of convergence and obtain more accurate solutions. We complement these optimized simultaneous sources and detectors by random simultaneous sources and detectors constrained to a complementary subspace. Our approach leads to solutions of the same quality as obtained using all sources and detectors but at a greatly reduced computational cost, as the number of large-scale linear systems to be solved is significantly reduced.

Toluwani Okunola, Mirjeta Pasha, Misha E. Kilmer

When reconstructing images from noisy measurements, such as in medical scans or scientific imaging, we face an inverse problem: recovering an unknown image from indirect, corrupted observations. These problems are typically ill-posed, meaning small amounts of noise can lead to inaccurate reconstructions. Regularization techniques address this by incorporating prior assumptions about the solution, such as smoothness or sparsity. However, standard methods often blur sharp edges--the boundaries between tissues or structures--losing critical detail. A powerful strategy for edge preservation is iterative reweighting, which solves a sequence of weighted subproblems with adaptively updated weights. Non-cumulative schemes derive weights from the current iterate alone and can be solved efficiently using the Recycled Majorization-Minimization Generalized Krylov Subspace method (RMM-GKS). The cumulative approach of Gazzola et al. progressively accumulates edge information across iterations, achieving superior edge preservation but at high computational cost. This work introduces CR-$\ell_q$-RMM-GKS, which combines cumulative edge detection with computational efficiency. We integrate Gazzola's cumulative weighting with RMM-GKS, which handles general $\ell_q$ penalties ($0 < q \le 2$), automatically selects regularization parameters, and recycles Krylov subspaces between iterations, reducing the nested structure to two levels. Numerical experiments in signal deblurring and tomography demonstrate that CR-$\ell_q$-RMM-GKS produces significantly sharper edge reconstructions than standard non-cumulative methods. In particular, CR-$\ell_1$-RMM-GKS outperforms both standard $\ell_1$ methods and CR-$\ell_2$-RMM-GKS, indicating that cumulative weighting and $\ell_1$ penalties are highly complementary.

Toluwani Okunola, Mirjeta Pasha, Misha E. Kilmer et al.

Many practical imaging systems suffer from uncertainty in acquisition geometry -- such as projection angles in computed tomography or sensor positions in photoacoustic tomography -- leading to nonlinear inverse problems that require joint estimation of both the image and the forward model parameters. Standard approaches that assume a known linear forward operator fail to account for these uncertainties, resulting in significant reconstruction artifacts. We propose a nonlinear recycled majorization-minimization generalized Krylov subspace (NL-RMM-GKS) framework for large-scale inverse problems with uncertain forward operators. The method extends MM-GKS to nonlinear settings by combining majorization-minimization for nonsmooth regularization with Krylov subspace projection and subspace recycling, ensuring bounded memory usage. Two complementary formulations are developed: an alternating minimization approach that alternates between image updates and Gauss-Newton parameter estimation, and a variable projection approach that eliminates the image variable and optimizes directly over the parameters using inexact inner solves. We further introduce streaming variants that process data sequentially, enabling reconstruction from large or dynamically acquired datasets without storing the full operator. For dynamic problems, we incorporate two temporal regularization strategies -- optical flow and anisotropic total variation -- as plug-in choices within the framework. We carry out rigorous numerical experiments in fan-beam computed tomography and photoacoustic tomography to demonstrate that our proposed framework achieves high-quality reconstructions with bounded memory requirements, making it suitable for large-scale dynamic imaging problems.

Osman Asif Malik, Shashanka Ubaru, Lior Horesh et al.

Many irregular domains such as social networks, financial transactions, neuron connections, and natural language constructs are represented using graph structures. In recent years, a variety of graph neural networks (GNNs) have been successfully applied for representation learning and prediction on such graphs. In many of the real-world applications, the underlying graph changes over time, however, most of the existing GNNs are inadequate for handling such dynamic graphs. In this paper we propose a novel technique for learning embeddings of dynamic graphs using a tensor algebra framework. Our method extends the popular graph convolutional network (GCN) for learning representations of dynamic graphs using the recently proposed tensor M-product technique. Theoretical results presented establish a connection between the proposed tensor approach and spectral convolution of tensors. The proposed method TM-GCN is consistent with the Message Passing Neural Network (MPNN) framework, accounting for both spatial and temporal message passing. Numerical experiments on real-world datasets demonstrate the performance of the proposed method for edge classification and link prediction tasks on dynamic graphs. We also consider an application related to the COVID-19 pandemic, and show how our method can be used for early detection of infected individuals from contact tracing data.

Elizabeth Newman, Misha E. Kilmer

In recent work (Soltani, Kilmer, Hansen, BIT 2016), an algorithm for non-negative tensor patch dictionary learning in the context of X-ray CT imaging and based on a tensor-tensor product called the $t$-product (Kilmer and Martin, 2011) was presented. Building on that work, in this paper, we use of non-negative tensor patch-based dictionaries trained on other data, such as facial image data, for the purposes of either compression or image deblurring. We begin with an analysis in which we address issues such as suitability of the tensor-based approach relative to a matrix-based approach, dictionary size and patch size to balance computational efficiency and qualitative representations. Next, we develop an algorithm that is capable of recovering non-negative tensor coefficients given a non-negative tensor dictionary. The algorithm is based on a variant of the Modified Residual Norm Steepest Descent method. We show how to augment the algorithm to enforce sparsity in the tensor coefficients, and note that the approach has broader applicability since it can be applied to the matrix case as well. We illustrate the surprising result that dictionaries trained on image data from one class can be successfully used to represent and compress image data from different classes and across different resolutions. Finally, we address the use of non-negative tensor dictionaries in image deblurring. We show that tensor treatment of the deblurring problem coupled with non-negative tensor patch dictionaries can give superior restorations as compared to standard treatment of the non-negativity constrained deblurring problem.

Sara Soltani, Misha E. Kilmer, Per Christian Hansen

We consider tomographic reconstruction using priors in the form of a dictionary learned from training images. The reconstruction has two stages: first we construct a tensor dictionary prior from our training data, and then we pose the reconstruction problem in terms of recovering the expansion coefficients in that dictionary. Our approach differs from past approaches in that a) we use a third-order tensor representation for our images and b) we recast the reconstruction problem using the tensor formulation. The dictionary learning problem is presented as a non-negative tensor factorization problem with sparsity constraints. The reconstruction problem is formulated in a convex optimization framework by looking for a solution with a sparse representation in the tensor dictionary. Numerical results show that our tensor formulation leads to very sparse representations of both the training images and the reconstructions due to the ability of representing repeated features compactly in the dictionary.

Oguz Semerci, Ning Hao, Misha E. Kilmer et al.

The development of energy selective, photon counting X-ray detectors allows for a wide range of new possibilities in the area of computed tomographic image formation. Under the assumption of perfect energy resolution, here we propose a tensor-based iterative algorithm that simultaneously reconstructs the X-ray attenuation distribution for each energy. We use a multi-linear image model rather than a more standard "stacked vector" representation in order to develop novel tensor-based regularizers. Specifically, we model the multi-spectral unknown as a 3-way tensor where the first two dimensions are space and the third dimension is energy. This approach allows for the design of tensor nuclear norm regularizers, which like its two dimensional counterpart, is a convex function of the multi-spectral unknown. The solution to the resulting convex optimization problem is obtained using an alternating direction method of multipliers (ADMM) approach. Simulation results shows that the generalized tensor nuclear norm can be used as a stand alone regularization technique for the energy selective (spectral) computed tomography (CT) problem and when combined with total variation regularization it enhances the regularization capabilities especially at low energy images where the effects of noise are most prominent.