In partial differential equations-based (PDE-based) inverse problems with many measurements, many... more In partial differential equations-based (PDE-based) inverse problems with many measurements, many largescale 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.
ABSTRACT Dissertation: G. W. Stewart III, “Some Topics in Numerical Analysis,” University of Tenn... more ABSTRACT Dissertation: G. W. Stewart III, “Some Topics in Numerical Analysis,” University of Tennessee. Published as Technical Report ORNL-4303, Oak Ridge National Laboratory, September 1968.
Parallel MR imaging is an effective approach to reduce MR image acquisition time. Non-uniform sub... more Parallel MR imaging is an effective approach to reduce MR image acquisition time. Non-uniform subsampling allows one to tailor the subsampling scheme for improved image quality at high acceleration factors. However, non-uniform subsampling precludes fast reconstruction schemes such as SENSE, and is more likely to require a regularized solution than reconstruction of uniformly subsampled data demands. This means that one needs to choose a good regularization parameter, typically requiring multiple expensive system solves. Here, we present an efficient LSQR-Hybrid algorithm which simultaneously addresses the need for rapid regularization parameter selection and fast reconstruction. This algorithm can reconstruct non-uniformly subsampled parallel MRI data, with automatic regularization and good image quality, in a time competitive with Cartesian SENSE.
In certain applications, the discretization of 2D integral equations can lead to system involving... more In certain applications, the discretization of 2D integral equations can lead to system involving matrices with block Toeplitz-Toeplitz block structure. Iterative Krylov subspace methods are sometimes employed to find regularized solutions to the related 2D discrete ill-posed problems; however, preconditioned which filter noise are needed to speed convergence to regularized solutions. We describe a preconditioning techniques based on the Toeplitz-type structure of the matrix which generalizes the approaches in (1) and (2) to take advantage of symmetry and real arithmetic operations. We use fast sine transforms to transform the original system to a system whose matrix has partially reconstructible Cauchy-like blocks. The preconditioner is a block diagonal, rank m approximation to this matrix, with Cauchy-like blocks each augmented by an identity of appropriate dimension. We note that the initialization cost is in general less than that for the similar 2D preconditioner in (2) which does not take advantage of symmetry. Several examples are given which illustrate the success of our preconditioned methods.
The imaging of shape perturbation and chromophore concentration using Diffuse Optical Tomography ... more The imaging of shape perturbation and chromophore concentration using Diffuse Optical Tomography (DOT) data can be mathematically described as an ill-posed and non-linear inverse problem. The reconstruction algorithm for hyperspectral data using a linearized Born model is prohibitively expensive, both in terms of computation and memory. We model the shape of the perturbation using parametric level-set approach (PaLS). We discuss novel computational strategies for reducing the computational cost based on a Krylov subspace approach for parameteric linear systems and a compression strategy for the parameter-to-observation map. We will demonstrate the validity of our approach by comparison with experiments.
We develop a parametric, shape-based image reconstruction algorithm for the joint reconstruction ... more We develop a parametric, shape-based image reconstruction algorithm for the joint reconstruction of the optical absorption and diffusion coefficients in the brain using diffuse optical tomographic data. Specifically, we study the recovery of the geometry of an unknown number of 2D closed contours located on a 2D manifold (the cortex) in 3-space. We describe an approach for a brain model in which we assume the existence of a one-to-one map from the surface of the cortex to a subset of the plane. We use a new, parametric level set approach to map shapes on the plane to structures on the cortex. Our optimization-based reconstruction algorithm evolves shapes on the plane while finding absorption and reduced scattering values inside each shape. Preliminary numerical simulation results show the promise of our approach.
Journal of Biomedical Materials Research Part A, Sep 24, 2012
Current statistical techniques for analyzing cellular alignment data in the fields of biomaterial... more Current statistical techniques for analyzing cellular alignment data in the fields of biomaterials and tissue engineering are limited due to heuristic and less quantitative approaches. For example, generally a cutoff degree limit (commonly 20 degrees), is arbitrarily defined within which cells are considered 'aligned.' The effectiveness of a patterned biomaterial in guiding the alignment of cells, such as neurons, is often critical to predict relationships between the biomaterial design and biological outcomes, both in vitro and in vivo. This becomes particularly important in the case of peripheral neurons, which require precise axon guidance to obtain successful regenerative outcomes. To address this issue, we have developed a protocol for processing cellular alignment data sets, which implicitly determines an 'angle of alignment.' This was accomplished as follows: cells 'aligning' with an underlying, anisotropic scaffold display uniformly distributed angles up to a cutoff point determined by how effective the biomaterial is in aligning cells. Therefore, this fact was then used to determine where an alignment angle data set diverges from a uniform distribution. This was accomplished by measuring the spacing between the collected, increasingly ordered angles and analyzing their underlying distributions using a normalized cumulative periodogram (NCP)-criterion. The proposed protocol offers a novel way to implicitly define cellular alignment, with respect to various anisotropic biomaterials. This method may also offer an alternative to assess cellular alignment, which could offer improved predictive measures related to biological outcomes. Further, the approach described can be used for a broad range of cell types grown on 2D surfaces, but would not be applicable to 3D scaffold systems in the present format.
In this paper, we develop a method for unsupervised clustering of twoway (matrix) data by combini... more In this paper, we develop a method for unsupervised clustering of twoway (matrix) data by combining two recent innovations from different fields: the Sparse Subspace Clustering (SSC) algorithm [10], which groups points coming from a union of subspaces into their respective subspaces, and the t-product [18], which was introduced to provide a matrix-like multiplication for third order tensors. Our algorithm is analogous to SSC in that an "affinity" between different data points is built using a sparse self-representation of the data. Unlike SSC, we employ the t-product in the self-representation. This allows us more flexibility in modeling; in fact, SSC is a special case of our method. When using the t-product, three-way arrays are treated as matrices whose elements (scalars) are n-tuples or tubes. Convolutions take the place of scalar multiplication. This framework allows us to embed the 2-D data into a vector-space-like structure called a free module over a commutative ring. These free modules retain many properties of complex inner-product spaces, and we leverage that to provide theoretical guarantees on our algorithm. We show that compared to vector-space counterparts, SSmC achieves higher accuracy and better able to cluster data with less preprocessing in some image clustering problems. In particular we show the performance of the proposed method on Weizmann face database, the Extended Yale B Face database and the MNIST handwritten digits database.
The discretization of integral equations can lead to systems involving symmetric Toeplitz matrice... more The discretization of integral equations can lead to systems involving symmetric Toeplitz matrices. We describe a preconditioning technique for the regularized solution of the related discrete ill-posed problem. We use discrete sine transforms to transform the system to one involving a Cauchy-like matrix. Based on the approachof Kilmer and O'Leary, the preconditioneris a symmetric, rank m approximation to the Cauchy-like matrix augmented by the identity. We shall show that if the kernel of the integral equation is smooth then the preconditioned matrix has two desirable properties; namely, the largest m magnitude eigenvalues are clustered around and bounded below by one, and that small magnitude eigenvalues remain small. We also show that the initialization cost is less than the initialization cost for the preconditioner introduced by Kilmer and O'Leary. Further, we describe a method for applying the preconditioner in O((n + 1) lg(n + 1)) operations when n+1 is a power of 2, and describe a variant of the MINRES algorithm to solve the symmetrically preconditioned problem. The preconditioned method is tested on two examples.
In partial differential equation-based (PDE-based) inverse problems with many measurements, many ... more In partial differential equation-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.
Loan were the organizers. About 35 people-computer scientists, mathematicians, and a broad range ... more Loan were the organizers. About 35 people-computer scientists, mathematicians, and a broad range of scholars who use tensor decompositions in their research-had come from eleven countries to participate in the weeklong workshop. Large group discussions and smaller break-out sessions were interleaved with invited talks, making for lively exchanges and a creative learning environment. Carla Martin, a graduate student at Cornell University, opened the week with an overview of multilinear algebra and tensor decompositions, beginning with the definition of a pth-order tensor A as a multiway array with p indices. A third-order tensor, for example, is written A = (a ijk) ∈ \ n 1× n 2× n 3. Workshop participants rated the workshop a success, admitting that they had raised at least as many questions as they answered.
Many ill-posed problems are solved using a discretization that results in a least squares problem... more Many ill-posed problems are solved using a discretization that results in a least squares problem or a linear system involving a Toeplitz matrix. The exact solution to such problems is often hopelessly contaminated by noise, since the discretized problem is quite ill-conditioned, and noise components in the approximate null-space dominate the solution vector. Therefore we seek an approximate solution that does not have large components in these directions. We use a preconditioned conjugate gradient algorithm to compute such a regularized solution. An orthogonal change of coordinates transforms the Toeplitz matrix to a Cauchy-like matrix, and we choose our preconditioner to be a low rank Cauchy-like matrix determined in the course of Gu's fast modi ed complete pivoting algorithm. We show that if the kernel of the ill-posed problem is smooth, then this preconditioner has desirable properties: the largest singular values of the preconditioned matrix are clustered around one, the smallest singular values, corresponding to the noise subspace, remain small, and the signal and noise spaces are relatively unmixed. The preconditioned algorithm costs only O(n lg n) operations per iteration for a problem with n variables. The e ectiveness of the preconditioner for ltering noise is demonstrated on three examples.
SIAM Journal on Matrix Analysis and Applications, 2000
In this paper we present three theorems which give insight into the regularizing properties of MI... more In this paper we present three theorems which give insight into the regularizing properties of MINRES. While our theory does not completely characterize the regularizing behavior of the algorithm, it provides a partial explanation of the observed behavior of the method. Unlike traditional attempts to explain the regularizing properties of Krylov subspace methods, our approach focuses on convergence properties of the residual rather than on convergence analysis of the harmonic Ritz values. The import of our analysis is illustrated by two examples. In particular, our theoretical and numerical results support the following important observation: in some circumstances the dimension of the optimal Krylov subspace can be much smaller than the number of the components of the truncated spectral solution that must be computed to attain comparable accuracy.
In this paper we present novel strategies for completion of 5D pre-stack seismic data, viewed as ... more In this paper we present novel strategies for completion of 5D pre-stack seismic data, viewed as a 5D tensor or as a set of 4D tensors across temporal frequencies. In contrast to existing complexity penalized algorithms for seismic data completion, which employ matrix analogues of tensor decompositions such as HOSVD or use overlapped Schatten norms from different unfoldings (or matricization) of the tensors, our approach uses a recently proposed decomposition called tensor SVD or tSVD for short, proposed in [Kilmer and Martin (2011)]. We show that seismic data exhibits low complexity under tSVD, i.e. is compressible under tSVD representation, and we subsequently propose a new complexity penalized algorithm for pre-stack seismic data completion under missing traces. This complexity measure which we call the Tensor Nuclear Norm (TNN) is motivated by algebraic properties of the tSVD. We test the performance of the proposed algorithms on synthetic and real data and show that missing data can be reliably recovered under heavy down-sampling.
In this paper, we introduce a method for multivariate function approximation using function evalu... more In this paper, we introduce a method for multivariate function approximation using function evaluations, Chebyshev polynomials, and tensor-based compression techniques via the Tucker format. We develop novel randomized techniques to accomplish the tensor compression, provide a detailed analysis of the computational costs, provide insight into the error of the resulting approximations, and discuss the benefits of the proposed approaches. We also apply the tensor-based function approximation to develop low-rank matrix approximations to kernel matrices that describe pairwise interactions between two sets of points; the resulting low-rank approximations are efficient to compute and store (the complexity is linear in the number of points). We have detailed numerical experiments on example problems involving multivariate function approximation, low-rank matrix approximations of kernel matrices involving well-separated clusters of sources and target points, and a global low-rank approximation of kernel matrices with an application to Gaussian processes. Keywords Multivariate function approximation • Tucker format • randomized algorithms • low-rank approximations • kernel methods A.K.
In partial differential equations-based (PDE-based) inverse problems with many measurements, many... more In partial differential equations-based (PDE-based) inverse problems with many measurements, many largescale 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.
ABSTRACT Dissertation: G. W. Stewart III, “Some Topics in Numerical Analysis,” University of Tenn... more ABSTRACT Dissertation: G. W. Stewart III, “Some Topics in Numerical Analysis,” University of Tennessee. Published as Technical Report ORNL-4303, Oak Ridge National Laboratory, September 1968.
Parallel MR imaging is an effective approach to reduce MR image acquisition time. Non-uniform sub... more Parallel MR imaging is an effective approach to reduce MR image acquisition time. Non-uniform subsampling allows one to tailor the subsampling scheme for improved image quality at high acceleration factors. However, non-uniform subsampling precludes fast reconstruction schemes such as SENSE, and is more likely to require a regularized solution than reconstruction of uniformly subsampled data demands. This means that one needs to choose a good regularization parameter, typically requiring multiple expensive system solves. Here, we present an efficient LSQR-Hybrid algorithm which simultaneously addresses the need for rapid regularization parameter selection and fast reconstruction. This algorithm can reconstruct non-uniformly subsampled parallel MRI data, with automatic regularization and good image quality, in a time competitive with Cartesian SENSE.
In certain applications, the discretization of 2D integral equations can lead to system involving... more In certain applications, the discretization of 2D integral equations can lead to system involving matrices with block Toeplitz-Toeplitz block structure. Iterative Krylov subspace methods are sometimes employed to find regularized solutions to the related 2D discrete ill-posed problems; however, preconditioned which filter noise are needed to speed convergence to regularized solutions. We describe a preconditioning techniques based on the Toeplitz-type structure of the matrix which generalizes the approaches in (1) and (2) to take advantage of symmetry and real arithmetic operations. We use fast sine transforms to transform the original system to a system whose matrix has partially reconstructible Cauchy-like blocks. The preconditioner is a block diagonal, rank m approximation to this matrix, with Cauchy-like blocks each augmented by an identity of appropriate dimension. We note that the initialization cost is in general less than that for the similar 2D preconditioner in (2) which does not take advantage of symmetry. Several examples are given which illustrate the success of our preconditioned methods.
The imaging of shape perturbation and chromophore concentration using Diffuse Optical Tomography ... more The imaging of shape perturbation and chromophore concentration using Diffuse Optical Tomography (DOT) data can be mathematically described as an ill-posed and non-linear inverse problem. The reconstruction algorithm for hyperspectral data using a linearized Born model is prohibitively expensive, both in terms of computation and memory. We model the shape of the perturbation using parametric level-set approach (PaLS). We discuss novel computational strategies for reducing the computational cost based on a Krylov subspace approach for parameteric linear systems and a compression strategy for the parameter-to-observation map. We will demonstrate the validity of our approach by comparison with experiments.
We develop a parametric, shape-based image reconstruction algorithm for the joint reconstruction ... more We develop a parametric, shape-based image reconstruction algorithm for the joint reconstruction of the optical absorption and diffusion coefficients in the brain using diffuse optical tomographic data. Specifically, we study the recovery of the geometry of an unknown number of 2D closed contours located on a 2D manifold (the cortex) in 3-space. We describe an approach for a brain model in which we assume the existence of a one-to-one map from the surface of the cortex to a subset of the plane. We use a new, parametric level set approach to map shapes on the plane to structures on the cortex. Our optimization-based reconstruction algorithm evolves shapes on the plane while finding absorption and reduced scattering values inside each shape. Preliminary numerical simulation results show the promise of our approach.
Journal of Biomedical Materials Research Part A, Sep 24, 2012
Current statistical techniques for analyzing cellular alignment data in the fields of biomaterial... more Current statistical techniques for analyzing cellular alignment data in the fields of biomaterials and tissue engineering are limited due to heuristic and less quantitative approaches. For example, generally a cutoff degree limit (commonly 20 degrees), is arbitrarily defined within which cells are considered 'aligned.' The effectiveness of a patterned biomaterial in guiding the alignment of cells, such as neurons, is often critical to predict relationships between the biomaterial design and biological outcomes, both in vitro and in vivo. This becomes particularly important in the case of peripheral neurons, which require precise axon guidance to obtain successful regenerative outcomes. To address this issue, we have developed a protocol for processing cellular alignment data sets, which implicitly determines an 'angle of alignment.' This was accomplished as follows: cells 'aligning' with an underlying, anisotropic scaffold display uniformly distributed angles up to a cutoff point determined by how effective the biomaterial is in aligning cells. Therefore, this fact was then used to determine where an alignment angle data set diverges from a uniform distribution. This was accomplished by measuring the spacing between the collected, increasingly ordered angles and analyzing their underlying distributions using a normalized cumulative periodogram (NCP)-criterion. The proposed protocol offers a novel way to implicitly define cellular alignment, with respect to various anisotropic biomaterials. This method may also offer an alternative to assess cellular alignment, which could offer improved predictive measures related to biological outcomes. Further, the approach described can be used for a broad range of cell types grown on 2D surfaces, but would not be applicable to 3D scaffold systems in the present format.
In this paper, we develop a method for unsupervised clustering of twoway (matrix) data by combini... more In this paper, we develop a method for unsupervised clustering of twoway (matrix) data by combining two recent innovations from different fields: the Sparse Subspace Clustering (SSC) algorithm [10], which groups points coming from a union of subspaces into their respective subspaces, and the t-product [18], which was introduced to provide a matrix-like multiplication for third order tensors. Our algorithm is analogous to SSC in that an "affinity" between different data points is built using a sparse self-representation of the data. Unlike SSC, we employ the t-product in the self-representation. This allows us more flexibility in modeling; in fact, SSC is a special case of our method. When using the t-product, three-way arrays are treated as matrices whose elements (scalars) are n-tuples or tubes. Convolutions take the place of scalar multiplication. This framework allows us to embed the 2-D data into a vector-space-like structure called a free module over a commutative ring. These free modules retain many properties of complex inner-product spaces, and we leverage that to provide theoretical guarantees on our algorithm. We show that compared to vector-space counterparts, SSmC achieves higher accuracy and better able to cluster data with less preprocessing in some image clustering problems. In particular we show the performance of the proposed method on Weizmann face database, the Extended Yale B Face database and the MNIST handwritten digits database.
The discretization of integral equations can lead to systems involving symmetric Toeplitz matrice... more The discretization of integral equations can lead to systems involving symmetric Toeplitz matrices. We describe a preconditioning technique for the regularized solution of the related discrete ill-posed problem. We use discrete sine transforms to transform the system to one involving a Cauchy-like matrix. Based on the approachof Kilmer and O'Leary, the preconditioneris a symmetric, rank m approximation to the Cauchy-like matrix augmented by the identity. We shall show that if the kernel of the integral equation is smooth then the preconditioned matrix has two desirable properties; namely, the largest m magnitude eigenvalues are clustered around and bounded below by one, and that small magnitude eigenvalues remain small. We also show that the initialization cost is less than the initialization cost for the preconditioner introduced by Kilmer and O'Leary. Further, we describe a method for applying the preconditioner in O((n + 1) lg(n + 1)) operations when n+1 is a power of 2, and describe a variant of the MINRES algorithm to solve the symmetrically preconditioned problem. The preconditioned method is tested on two examples.
In partial differential equation-based (PDE-based) inverse problems with many measurements, many ... more In partial differential equation-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.
Loan were the organizers. About 35 people-computer scientists, mathematicians, and a broad range ... more Loan were the organizers. About 35 people-computer scientists, mathematicians, and a broad range of scholars who use tensor decompositions in their research-had come from eleven countries to participate in the weeklong workshop. Large group discussions and smaller break-out sessions were interleaved with invited talks, making for lively exchanges and a creative learning environment. Carla Martin, a graduate student at Cornell University, opened the week with an overview of multilinear algebra and tensor decompositions, beginning with the definition of a pth-order tensor A as a multiway array with p indices. A third-order tensor, for example, is written A = (a ijk) ∈ \ n 1× n 2× n 3. Workshop participants rated the workshop a success, admitting that they had raised at least as many questions as they answered.
Many ill-posed problems are solved using a discretization that results in a least squares problem... more Many ill-posed problems are solved using a discretization that results in a least squares problem or a linear system involving a Toeplitz matrix. The exact solution to such problems is often hopelessly contaminated by noise, since the discretized problem is quite ill-conditioned, and noise components in the approximate null-space dominate the solution vector. Therefore we seek an approximate solution that does not have large components in these directions. We use a preconditioned conjugate gradient algorithm to compute such a regularized solution. An orthogonal change of coordinates transforms the Toeplitz matrix to a Cauchy-like matrix, and we choose our preconditioner to be a low rank Cauchy-like matrix determined in the course of Gu's fast modi ed complete pivoting algorithm. We show that if the kernel of the ill-posed problem is smooth, then this preconditioner has desirable properties: the largest singular values of the preconditioned matrix are clustered around one, the smallest singular values, corresponding to the noise subspace, remain small, and the signal and noise spaces are relatively unmixed. The preconditioned algorithm costs only O(n lg n) operations per iteration for a problem with n variables. The e ectiveness of the preconditioner for ltering noise is demonstrated on three examples.
SIAM Journal on Matrix Analysis and Applications, 2000
In this paper we present three theorems which give insight into the regularizing properties of MI... more In this paper we present three theorems which give insight into the regularizing properties of MINRES. While our theory does not completely characterize the regularizing behavior of the algorithm, it provides a partial explanation of the observed behavior of the method. Unlike traditional attempts to explain the regularizing properties of Krylov subspace methods, our approach focuses on convergence properties of the residual rather than on convergence analysis of the harmonic Ritz values. The import of our analysis is illustrated by two examples. In particular, our theoretical and numerical results support the following important observation: in some circumstances the dimension of the optimal Krylov subspace can be much smaller than the number of the components of the truncated spectral solution that must be computed to attain comparable accuracy.
In this paper we present novel strategies for completion of 5D pre-stack seismic data, viewed as ... more In this paper we present novel strategies for completion of 5D pre-stack seismic data, viewed as a 5D tensor or as a set of 4D tensors across temporal frequencies. In contrast to existing complexity penalized algorithms for seismic data completion, which employ matrix analogues of tensor decompositions such as HOSVD or use overlapped Schatten norms from different unfoldings (or matricization) of the tensors, our approach uses a recently proposed decomposition called tensor SVD or tSVD for short, proposed in [Kilmer and Martin (2011)]. We show that seismic data exhibits low complexity under tSVD, i.e. is compressible under tSVD representation, and we subsequently propose a new complexity penalized algorithm for pre-stack seismic data completion under missing traces. This complexity measure which we call the Tensor Nuclear Norm (TNN) is motivated by algebraic properties of the tSVD. We test the performance of the proposed algorithms on synthetic and real data and show that missing data can be reliably recovered under heavy down-sampling.
In this paper, we introduce a method for multivariate function approximation using function evalu... more In this paper, we introduce a method for multivariate function approximation using function evaluations, Chebyshev polynomials, and tensor-based compression techniques via the Tucker format. We develop novel randomized techniques to accomplish the tensor compression, provide a detailed analysis of the computational costs, provide insight into the error of the resulting approximations, and discuss the benefits of the proposed approaches. We also apply the tensor-based function approximation to develop low-rank matrix approximations to kernel matrices that describe pairwise interactions between two sets of points; the resulting low-rank approximations are efficient to compute and store (the complexity is linear in the number of points). We have detailed numerical experiments on example problems involving multivariate function approximation, low-rank matrix approximations of kernel matrices involving well-separated clusters of sources and target points, and a global low-rank approximation of kernel matrices with an application to Gaussian processes. Keywords Multivariate function approximation • Tucker format • randomized algorithms • low-rank approximations • kernel methods A.K.
Uploads
Papers by Misha Kilmer