Analysis of eigenvalue condition numbers for a class of randomized numerical methods for singular matrix pencils.

The numerical solution of the generalized eigenvalue problem for a singular matrix pencil is challenging due to the discontinuity of its eigenvalues. Classically, such problems are addressed by first extracting the regular part through the staircase form and then applying a standard solver, such as the QZ algorithm, to that regular part. Recently, several novel approaches have been proposed to transform the singular pencil into a regular pencil by relatively simple randomized modifications. In this work, we analyze three such methods by Hochstenbach, Mehl, and Plestenjak that modify, project, or augment the pencil using random matrices. All three methods rely on the normal rank and do not alter the finite eigenvalues of the original pencil. We show that the eigenvalue condition numbers of the transformed pencils are unlikely to be much larger than the δ -weak eigenvalue condition numbers, introduced by Lotz and Noferini, of the original pencil. This not only indicates favorable numerical stability but also reconfirms that these condition numbers are a reliable criterion for detecting simple finite eigenvalues. We also provide evidence that, from a numerical stability perspective, the use of complex instead of real random matrices is preferable even for real singular matrix pencils and real eigenvalues. As a side result, we provide sharp left tail bounds for a product of two independent random variables distributed with the generalized beta distribution of the first kind or Kumaraswamy distribution.

Singular-Value-Decomposition-Based Matrix Surgery.

This paper is motivated by the need to stabilise the impact of deep learning (DL) training for medical image analysis on the conditioning of convolution filters in relation to model overfitting and robustness. We present a simple strategy to reduce square matrix condition numbers and investigate its effect on the spatial distributions of point clouds of well- and ill-conditioned matrices. For a square matrix, the SVD surgery strategy works by: (1) computing its singular value decomposition (SVD), (2) changing a few of the smaller singular values relative to the largest one, and (3) reconstructing the matrix by reverse SVD. Applying SVD surgery on CNN convolution filters during training acts as spectral regularisation of the DL model without requiring the learning of extra parameters. The fact that the further away a matrix is from the non-invertible matrices, the higher its condition number is suggests that the spatial distributions of square matrices and those of their inverses are correlated to their condition number distributions. We shall examine this assertion empirically by showing that applying various versions of SVD surgery on point clouds of matrices leads to bringing their persistent diagrams (PDs) closer to the matrices of the point clouds of their inverses.

Ultra-limited-angle CT image reconstruction algorithm based on reweighting and edge-preserving.

BACKGROUND: Ultra-limited-angle image reconstruction problem with a limited-angle scanning range less than or equal to π2 is severely ill-posed. Due to the considerably large condition number of a linear system for image reconstruction, it is extremely challenging to generate a valid reconstructed image by traditional iterative reconstruction algorithms. OBJECTIVE: To develop and test a valid ultra-limited-angle CT image reconstruction algorithm. METHODS: We propose a new optimized reconstruction model and Reweighted Alternating Edge-preserving Diffusion and Smoothing algorithm in which a reweighted method of improving the condition number is incorporated into the idea of AEDS image reconstruction algorithm. The AEDS algorithm utilizes the property of image sparsity to improve partially the results. In experiments, the different algorithms (the Pre-Landweber, AEDS algorithms and our algorithm) are used to reconstruct the Shepp-Logan phantom from the simulated projection data with noises and the flat object with a large ratio between length and width from the real projection data. PSNR and SSIM are used as the quantitative indices to evaluate quality of reconstructed images. RESULTS: Experiment results showed that for simulated projection data, our algorithm improves PSNR and SSIM from 22.46db to 39.38db and from 0.71 to 0.96, respectively. For real projection data, our algorithm yields the highest PSNR and SSIM of 30.89db and 0.88, which obtains a valid reconstructed result. CONCLUSIONS: Our algorithm successfully combines the merits of several image processing and reconstruction algorithms. Thus, our new algorithm outperforms significantly other two algorithms and is valid for ultra-limited-angle CT image reconstruction.

A Critical Evaluation and Modification of the Padé-Laplace Method for Deconvolution of Viscoelastic Spectra.

This paper shows that using the Padé-Laplace (PL) method for deconvolution of multi-exponential functions (stress relaxation of polymers) can produce ill-conditioned systems of equations. Analysis of different sets of generated data points from known multi-exponential functions indicates that by increasing the level of Padé approximants, the condition number of a matrix whose entries are coefficients of a Taylor series in the Laplace space grows rapidly. When higher levels of Padé approximants need to be computed to achieve stable modes for separation of exponentials, the problem of generating matrices with large condition numbers becomes more pronounced. The analysis in this paper discusses the origin of ill-posedness of the PL method and it was shown that ill-posedness may be regularized by reconstructing the system of equations and using singular value decomposition (SVD) for computation of the Padé table. Moreover, it is shown that after regularization, the PL method can deconvolute the exponential decays even when the input parameter of the method is out of its optimal range.

A preconditioned landweber iteration scheme for the limited-angle image reconstruction.

BACKGROUND: The limited-angle reconstruction problem is of both theoretical and practical importance. Due to the severe ill-posedness of the problem, it is very challenging to get a valid reconstructed result from the known small limited-angle projection data. The theoretical ill-posedness leads the normal equation AT Ax = AT b of the linear system derived by discretizing the Radon transform to be severely ill-posed, which is quantified as the large condition number of AT A. OBJECTIVE: To develop and test a new valid algorithm for improving the limited-angle image reconstruction with the known appropriately small angle range from [0,π3]∼[0,π2]. METHODS: We propose a reweighted method of improving the condition number of AT Ax = AT b and the corresponding preconditioned Landweber iteration scheme. The weight means multiplying AT Ax = AT b by a matrix related to AT A, and the weighting process is repeated multiple times. In the experiment, the condition number of the coefficient matrix in the reweighted linear system decreases monotonically to 1 as the weighting times approaches infinity. RESULTS: The numerical experiments showed that the proposed algorithm is significantly superior to other iterative algorithms (Landweber, Cimmino, NWL-a and AEDS) and can reconstruct a valid image from the known appropriately small angle range. CONCLUSIONS: The proposed algorithm is effective for the limited-angle reconstruction problem with the known appropriately small angle range.

Efficient newton-raphson/singular value decomposition-based optimization scheme with dynamically updated critical condition number for rapid convergence of weighted histogram analysis method equations.

Selection of the successful optimization strategy is an essential part of solving numerous practical problems yet often is a nontrivial task, especially when a function to be optimized is multidimensional and involves statistical data. Here we propose a robust optimization scheme, referred to as NR/SVD-Cdyn, which is based on a combination of the Newton-Raphson (NR) method along with singular value decomposition (SVD), and demonstrate its performance by numerically solving a system of the weighted histogram analysis method equations. Our results show significant improvement over the direct iteration and conventional NR optimization methods. The proposed scheme is universal and could be used for solving various optimization problems in the field of computational chemistry such as parameter fitting for the methods of molecular mechanics and semiempirical quantum-mechanical methods. © 2019 Wiley Periodicals, Inc.

Scalable estimation and regularization for the logistic normal multinomial model.

Clustered multinomial data are prevalent in a variety of applications such as microbiome studies, where metagenomic sequencing data are summarized as multinomial counts for a large number of bacterial taxa per subject. Count normalization with ad hoc zero adjustment tends to result in poor estimates of abundances for taxa with zero or small counts. To account for heterogeneity and overdispersion in such data, we suggest using the logistic normal multinomial (LNM) model with an arbitrary correlation structure to simultaneously estimate the taxa compositions by borrowing information across subjects. We overcome the computational difficulties in high dimensions by developing a stochastic approximation EM algorithm with Hamiltonian Monte Carlo sampling for scalable parameter estimation in the LNM model. The ill-conditioning problem due to unstructured covariance is further mitigated by a covariance-regularized estimator with a condition number constraint. The advantages of the proposed methods are illustrated through simulations and an application to human gut microbiome data.

A multicontrast imaging method using steady-state free precession with alternating RF flip angles.

PURPOSE: To obtain multicontrast images including fat-suppressed contrast image, a novel multicontrast imaging method using an SSFP sequence with alternating RF flip angles is proposed. METHODS: The proposed method uses the balanced SSFP sequence with 2 flip angles. In general, the conventional balanced SSFP sequence has its own unique contrast, which combines both FID signal and echo signal under a steady-state condition. By using alternating RF flip angles and RF phase cycling, various image contrasts weighted by proton density, T1 , and T2 can be obtained. The proposed method offers multicontrast images with fat suppression by using the combination of 2 images obtained just after alternating RF pulses, respectively. RESULTS: As demonstrated by simulations, phantom and in vivo experiments, the proposed method provides multicontrast knee images including fat-suppressed contrast images. CONCLUSION: The proposed method can be a useful tool for clinical diagnosis, such as the cartilage segmentation and the fast screening of lesions.

Optimal Shrinkage of Eigenvalues in the Spiked Covariance Model.

We show that in a common high-dimensional covariance model, the choice of loss function has a profound effect on optimal estimation. In an asymptotic framework based on the Spiked Covariance model and use of orthogonally invariant estimators, we show that optimal estimation of the population covariance matrix boils down to design of an optimal shrinker η that acts elementwise on the sample eigenvalues. Indeed, to each loss function there corresponds a unique admissible eigenvalue shrinker η* dominating all other shrinkers. The shape of the optimal shrinker is determined by the choice of loss function and, crucially, by inconsistency of both eigenvalues and eigenvectors of the sample covariance matrix. Details of these phenomena and closed form formulas for the optimal eigenvalue shrinkers are worked out for a menagerie of 26 loss functions for covariance estimation found in the literature, including the Stein, Entropy, Divergence, Fréchet, Bhattacharya/Matusita, Frobenius Norm, Operator Norm, Nuclear Norm and Condition Number losses.

A novel measure of reliability in Diffusion Tensor Imaging after data rejections due to subject motion.

Diffusion Tensor Imaging (DTI) is commonly challenged by subject motion during data acquisition, which often leads to corrupted image data. Currently used procedure in DTI analysis is to correct or completely reject such data before tensor estimations, however assessing the reliability and accuracy of the estimated tensor in such situations has evaded previous studies. This work aims to define the loss of data accuracy with increasing image rejections, and to define a robust method for assessing reliability of the result at voxel level. We carried out simulations of every possible sub-scheme (N=1,073,567,387) of Jones30 gradient scheme, followed by confirming the idea with MRI data from four newborn and three adult subjects. We assessed the relative error of the most commonly used tensor estimates for DTI and tractography studies, fractional anisotropy (FA) and the major orientation vector (V1), respectively. The error was estimated using two measures, the widely used electric potential (EP) criteria as well as the rotationally variant condition number (CN). Our results show that CN and EP are comparable in situations with very few rejections, but CN becomes clearly more sensitive to depicting errors when more gradient vectors and images were rejected. The error in FA and V1 was also found depend on the actual FA level in the given voxel; low actual FA levels were related to high relative errors in the FA and V1 estimates. Finally, the results were confirmed with clinical MRI data. This showed that the errors after rejections are, indeed, inhomogeneous across brain regions. The FA and V1 errors become progressively larger when moving from the thick white matter bundles towards more superficial subcortical structures. Our findings suggest that i) CN is a useful estimator of data reliability at voxel level, and ii) DTI preprocessing with data rejections leads to major challenges when assessing brain tissue with lower FA levels, such as all newborn brain, as well as the adult superficial, subcortical areas commonly traced in precise connectivity analyses between cortical regions.

Signal Conditioning for the Kalman Filter: Application to Satellite Attitude Estimation with Magnetometer and Sun Sensors.

Most satellites use an on-board attitude estimation system, based on available sensors. In the case of low-cost satellites, which are of increasing interest, it is usual to use magnetometers and Sun sensors. A Kalman filter is commonly recommended for the estimation, to simultaneously exploit the information from sensors and from a mathematical model of the satellite motion. It would be also convenient to adhere to a quaternion representation. This article focuses on some problems linked to this context. The state of the system should be represented in observable form. Singularities due to alignment of measured vectors cause estimation problems. Accommodation of the Kalman filter originates convergence difficulties. The article includes a new proposal that solves these problems, not needing changes in the Kalman filter algorithm. In addition, the article includes assessment of different errors, initialization values for the Kalman filter; and considers the influence of the magnetic dipole moment perturbation, showing how to handle it as part of the Kalman filter framework.

Paradoxical effect of the signal-to-noise ratio of GRAPPA calibration lines: A quantitative study.

PURPOSE: Intuitively, GRAPPA auto-calibration signal (ACS) lines with higher signal-to-noise ratio (SNR) may be expected to boost the accuracy of kernel estimation and increase the SNR of GRAPPA reconstructed images. Paradoxically, Sodickson and his colleagues pointed out that using ACS lines with high SNR may actually lead to lower SNR in the GRAPPA reconstructed images. A quantitative study of how the noise in the ACS lines affects the SNR of the GRAPPA reconstructed images is presented. METHODS: In a phantom, the singular values of the GRAPPA encoding matrix and the root-mean-square error of GRAPPA reconstruction were evaluated using multiple sets of ACS lines with variant SNR. In volunteers, ACS lines with high and low SNR were estimated, and the SNR of corresponding TGRAPPA reconstructed images was evaluated. RESULTS: We show that the condition number of the GRAPPA kernel estimation equations is proportional to the SNR of the ACS lines. In dynamic image series reconstructed with TGRAPPA, high SNR ACS lines result in reduced SNR if appropriate regularization is not applied. CONCLUSION: Noise has a similar effect to Tikhonov regularization. Without appropriate regularization, a GRAPPA kernel estimated from ACS lines with higher SNR amplifies random noise in the GRAPPA reconstruction. Magn Reson Med 74:231-239, 2015. © 2014 Wiley Periodicals, Inc.

System matrix analysis for sparse-view iterative image reconstruction in X-ray CT.

Iterative image reconstruction (IIR) with sparsity-exploiting methods, such as total variation (TV) minimization, used for investigations in compressive sensing (CS) claim potentially large reductions in sampling requirements. Quantifying this claim for computed tomography (CT) is non-trivial, as both the singularity of undersampled reconstruction and the sufficient view number for sparse-view reconstruction are ill-defined. In this paper, the singular value decomposition method is used to study the condition number and singularity of the system matrix and the regularized matrix. An estimation method of the empirical lower bound is proposed, which is helpful for estimating the number of projection views required for exact reconstruction. Simulation studies show that the singularity of the system matrices for different projection views is effectively reduced by regularization. Computing the condition number of a regularized matrix is necessary to provide a reference for evaluating the singularity and recovery potential of reconstruction algorithms using regularization. The empirical lower bound is helpful for estimating the projections view number with a sparse reconstruction algorithm.

Stable Estimation of a Covariance Matrix Guided by Nuclear Norm Penalties.

Estimation of a covariance matrix or its inverse plays a central role in many statistical methods. For these methods to work reliably, estimated matrices must not only be invertible but also well-conditioned. The current paper introduces a novel prior to ensure a well-conditioned maximum a posteriori (MAP) covariance estimate. The prior shrinks the sample covariance estimator towards a stable target and leads to a MAP estimator that is consistent and asymptotically efficient. Thus, the MAP estimator gracefully transitions towards the sample covariance matrix as the number of samples grows relative to the number of covariates. The utility of the MAP estimator is demonstrated in two standard applications - discriminant analysis and EM clustering - in this sampling regime.

Immunity to Increasing Condition Numbers of Linear Superiorization versus Linear Programming.

Given a family of linear constraints and a linear objective function one can consider whether to apply a Linear Programming (LP) algorithm or use a Linear Superiorization (LinSup) algorithm on this data. In the LP methodology one aims at finding a point that fulfills the constraints and has the minimal value of the objective function over these constraints. The Linear Superiorization approach considers the same data as linear programming problems but instead of attempting to solve those with linear optimization methods it employs perturbation resilient feasibility-seeking algorithms and steers them toward feasible points with reduced (not necessarily minimal) objective function values. Previous studies compared LP and LinSup in terms of their respective outputs and the resources they use. In this paper we investigate these two approaches in terms of their sensitivity to condition numbers of the system of linear constraints. Condition numbers are a measure for the impact of deviations in the input data on the output of a problem and, in particular, they describe the factor of error propagation when given wrong or erroneous data. Therefore, the ability of LP and LinSup to cope with increased condition numbers, thus with ill-posed problems, is an important matter to consider which was not studied until now. We investigate experimentally the advantages and disadvantages of both LP and LinSup on examplary problems of linear programming with multiple condition numbers and different problem dimensions.

An improved composite ship magnetic field model with ellipsoid and magnetic dipole arrays.

In order to simultaneously maintain the ship magnetic field modeling accuracy, reduce the number of coefficient matrix conditions and the model computational complexity, an improved composite model is designed by introducing the magnetic dipole array model with a single-axis magnetic moment on the basis of the hybrid ellipsoid and magnetic dipole array model. First, the improved composite model of the ship's magnetic field is established based on the magnetic dipole array model with 3-axis magnetic moment, the magnetic dipole array model with only x-axis magnetic moment, and the ellipsoid model. Secondly, the set of equations for calculating the magnetic moments of the composite model is established, and for the problem of solving the pathological set of equations, the least-squares estimation, stepwise regression method, Tikhonov, and truncated singular value decomposition regularization methods are introduced in terms of the magnetic field, and generalized cross-validation is used to solve the optimal regularization parameters. Finally, a ship model test is designed to compare and analyze the effectiveness of the composite and hybrid models in four aspects: the number of coefficient matrix conditions of the model equation set, the relative error of magnetic field fitting, the relative error of magnetic field extrapolation, and the computational time complexity. The modeling results based on the ship model test data show that the composite model can be used for modeling the magnetic field of ships, and compared with the hybrid model, it reduces the number of coefficient matrix conditions and improves the computational efficiency on the basis of retaining a higher modeling accuracy, and it can be effectively applied in related scientific research and engineering.

Adaptive Fitting of Linear Mixed-Effects Models with Correlated Random-effects.

Linear mixed-effects model has been widely used in longitudinal data analyses. In practice, the fitting algorithm can fail to converge due to boundary issues of the estimated random-effects covariance matrix G, i.e., being near-singular, non-positive definite, or both. Current available algorithms are not computationally optimal because the condition number of matrix G is unnecessarily increased when the random-effects correlation estimate is not zero. We propose an adaptive fitting (AF) algorithm using an optimal linear transformation of the random-effects design matrix. It is a data-driven adaptive procedure, aiming at reducing subsequent random-effects correlation estimates down to zero in the optimal transformed estimation space. Simulations show that AF significantly improves the convergent properties, especially under small sample size, relative large noise and high correlation settings. One real data for Insulin-like Growth Factor (IGF) protein is used to illustrate the application of this algorithm implemented with software package R (nlme).

Band selection for mapping chromophores of skin tissue.

A numerical approach has been proposed to identify bands for optimally estimating the concentration of three types of viable chromophores within biological tissue. The bands are determined according to the condition number of absorption matrix associated with the attenuation coefficients of chromophores. The effectiveness of different sets of selected band combination was verified by using the spectral reflectance images of skin tissue acquired from standard forearm vascular occlusion tests via a spectroradiometer. Experimental results demonstrated that the concentration of chromophores within skin tissue could be estimated correctly and robustly only when the bands were deliberately selected.

Experimental demonstration of improved magnetorelaxometry imaging performance using optimized coil configurations.

BACKGROUND: Magnetorelaxometry imaging is an experimental imaging technique capable of reconstructing magnetic nanoparticle distributions inside a volume noninvasively and with high specificity. Thus, magnetorelaxometry imaging is a promising candidate for monitoring a number of therapeutical approaches that employ magnetic nanoparticles, such as magnetic drug targeting and magnetic hyperthermia, to guarantee their safety and efficacy. Prior to a potential clinical application of this imaging modality, it is necessary to optimize magnetorelaxometry imaging systems to produce reliable imaging results and to maximize the reconstruction accuracy of the magnetic nanoparticle distributions. Multiple optimization approaches were already applied throughout a number of simulation studies, all of which yielded increased imaging qualities compared to intuitively designed measurement setups. PURPOSE: None of these simulative approaches was conducted in practice such that it still remains unclear if the theoretical results are achievable in an experimental setting. In this study, we demonstrate the technical feasibility and the increased reconstruction accuracy of optimized coil configurations in two distinct magnetorelaxometry setups. METHODS: The electromagnetic coil positions and radii of a cuboidal as well as a cylindrical magnetorelaxometry imaging setup are optimized by minimizing the system matrix condition numbers of their corresponding linear forward models. The optimized coil configurations are manufactured alongside with two regular coil grids. Magnetorelaxometry measurements of three cuboidal and four cylindrical magnetic nanoparticle phantoms are conducted, and the resulting reconstruction qualities of the optimized and the regular coil configurations are compared. RESULTS: The computed condition numbers of the optimized coil configurations are approximately one order of magnitude lower compared to the regular coil grids. The reconstruction results show that for both setups, every phantom is recovered more accurately by the optimized coil configurations compared to the regular coil grids. Additionally, the optimized coil configurations yield better signal qualities. CONCLUSIONS: The presented experimental study provides a proof of the practicality and the efficacy of optimizing magnetorelaxometry imaging systems with respect to the condition numbers of their system matrices, previously only demonstrated in simulations. From the promising results of our study, we infer that the minimization of the system matrix condition number will also enable the practical optimization of other design parameters of magnetorelaxometry imaging setups (e.g., sensor configuration, coil currents, etc.) in order to improve the achievable reconstruction qualities even further, eventually paving the way towards clinical application of this imaging modality.

Modeling the impact of point mutations on the stability of proteins.

A new method has been introduced which allows us to determine the stability of protein complexes with point changes of amino acid residues that also take into account the three-dimensional structure of the complex. This formulated and proven theorem is aimed at determining the criterion for the stability of protein molecules. The algorithm and software package were developed for analyzing protein interactions, taking into account their three-dimensional structure from the PDB database.