Small sample statistical condition estimation for the total least squares problem

In this paper, under the genericity condition, we study the condition estimation of the total least squares (TLS) problem based on small sample condition estimation (SCE), which can be incorporated into the direct solver for the TLS problem via the singular value decomposition (SVD) of the augmented matrix [A, b]. Our proposed condition estimation algorithms are efficient for the small and medium size TLS problem because they utilize the computed SVD of [A, b] during the numerical solution to the TLS problem. Numerical examples illustrate the reliability of the algorithms. Both normwise and componentwise perturbations are considered. Moreover, structured condition estimations are investigated for the structured TLS problem.     

Accurate SVDs of weakly diagonally dominant M-matrices

Summary. We present a new O(n³) algorithm which computes the SVD of a weakly diagonally dominant M-matrix to high relative accuracy. The algorithm takes as an input the offdiagonal entries of the matrix and its row sums.Mathematics Subject Classification (1991): 65F15Revised version received September 19, 2003This material is based in part upon work supported by the LLNL Memorandum Agreement No. B504962 under DOE Contract No. W-7405-ENG-48, DOE Grants No. DE-FG03-94ER25219, DE-FC03-98ER25351 and DE-FC02-01ER25478, NSF Grant No. ASC-9813362, and Cooperative Agreement No. ACI-9619020.     

Optimal recovery of 3D X-ray tomographic data via shearlet decomposition

This paper introduces a new decomposition of the 3D X-ray transform based on the shearlet representation, a multiscale directional representation which is optimally efficient in handling 3D data containing edge singularities. Using this decomposition, we derive a highly effective reconstruction algorithm yielding a near-optimal rate of convergence in estimating piecewise smooth objects from 3D X-ray tomographic data which are corrupted by white Gaussian noise. This algorithm is achieved by applying a thresholding scheme on the 3D shearlet transform coefficients of the noisy data which, for a given noise level ε, can be tuned so that the estimator attains the essentially optimal mean square error rate O(log(ε ^???1)ε ^2/3), as ε→0. This is the first published result to achieve this type of error estimate, outperforming methods based on Wavelet-Vaguelettes decomposition and on SVD, which can only achieve MSE rates of O(ε ^1/2) and O(ε ^1/3), respectively.     

Simultaneous singular value decomposition

We consider the following problem: Given a set of m×n real (or complex) matrices A₁,…,A_N, find an m×m orthogonal (or unitary) matrix P and an n×n orthogonal (or unitary) matrix Q such that P^*A₁Q,…,P^*A_NQ are in a common block-diagonal form with possibly rectangular diagonal blocks. We call this the simultaneous singular value decomposition (simultaneous SVD). The name is motivated by the fact that the special case with N=1, where a single matrix is given, reduces to the ordinary SVD. With the aid of the theory of *-algebra and bimodule it is shown that a finest simultaneous SVD is uniquely determined. An algorithm is proposed for finding the finest simultaneous SVD on the basis of recent algorithms of Murota-Kanno-Kojima-Kojima and Maehara-Murota for simultaneous block-diagonalization of square matrices under orthogonal (or unitary) similarity.     

Simulation of the interaction of light with 3-D metallic nanostructures using a proper orthogonal decomposition-Galerkin reduced-order discontinuous Galerkin time-domain method

In this artice, we report on a reduced-order model (ROM) based on the proper orthogonal decomposition (POD) technique for the system of 3-D time-domain Maxwell's equations coupled to a Drude dispersion model, which is employed to describe the interaction of light with nanometer scale metallic structures. By using the singular value decomposition (SVD) method, the POD basis vectors are extracted offline from the snapshots produced by a high order discontinuous Galerkin time-domain (DGTD) solver. With a Galerkin projection and a second order leap-frog (LF₂) time discretization, a discrete ROM is constructed. The stability condition of the ROM is then analyzed. In particular, when the boundary is a perfect electric conductor condition, the global energy of the ROM is bounded, which is consistent with the characteristics of global energy in the DGTD method. It is shown that the ROM based on Galerkin projection can maintain the stability characteristics of the original high dimensional model. Numerical experiments are presented to verify the accuracy, demonstrate the capabilities of the POD-based ROM and assess its efficiency for 3-D nanophotonic problems.     

Computational uncertainty and the application of a high-performance multiple precision scheme to obtaining the correct reference solution of Lorenz equations

The computational uncertainty principle (CUP) is applied to explain the experimental formulae of the critical time of decoupling for Lorenz equations (LEs). We apply the multiple precision (MP) library in obtaining the long-time solution of LEs, and based on the classic Taylor scheme, we developed a high-performance parallel Taylor solver to do the computation. The new solver is several hundreds times faster than the reported solvers developed in MATHEMATICA software, and it has the ability to yield longer solutions of LEs, up to t ∼ 10⁴ LTU (Lorenz time unit). Further, we notice that the two computation processes with different precisions or orders will produce the reliable correct reference solutions before they have a significant difference. According to this property we propose an approach for maintaining the correct numerical solution. The new solver and the solution validation approach are used to identify and correct an erroneous solution reported in a previous study.     

Performance comparison of the NWF and DC methods for implementing High-Resolution schemes in a fully coupled incompressible flow solver

This paper reports on the use of the Normalized Weighting Factor (NWF) method and the Deferred Correction (DC) approach for the implementation of High Resolution (HR) convective schemes in an implicit, fully coupled, pressure-based flow solver. Four HR schemes are realized within the framework of the NWF and DC methods and employed to solve the following three laminar flow problems: (i) lid-driven flow in a square cavity, (ii) sudden expansion in a square cavity, and (iii) flow in a planar T-junction, over three grid systems with sizes of 10⁴, 5 × 10⁴, and 3 × 10⁵ control volumes. The merit of both approaches is demonstrated by comparing the computational costs required to solve these problems using the various HR schemes on the different grid systems. Whereas previous attempts to use the NWF method in a segregated flow solver failed to produce converged solutions, current results clearly demonstrate that both methods are suitable for utilization in a coupled flow solver. In terms of CPU efficiency, there is no global and consistent superiority of any method over another even though the DC method outperformed the NWF method in two of the three test problems solved.     

Statistical theory of shape under elliptical models and singular value decompositions

The size-and-shape and shape distributions based on non-central and non-isotropic elliptical distributions are derived in this paper by using the singular value decomposition (SVD). The general densities require the computation of new integrals involving zonal polynomials. The invariance of the central shape distribution is also proved. Finally, some particular densities are applied in a classical data of Biology, and the inference based on exact distributions is performed after choosing the best model by using a modified BIC^∗ criterion.     

Solution of large‐scale weighted least‐squares problems

A sequence of least‐squares problems of the form min_y∥G^1/2(A^T y?h)∥₂, where G is an n×n positive‐definite diagonal weight matrix, and A an m×n (m?n) sparse matrix with some dense columns; has many applications in linear programming, electrical networks, elliptic boundary value problems, and structural analysis. We suggest low‐rank correction preconditioners for such problems, and a mixed solver (a combination of a direct solver and an iterative solver). The numerical results show that our technique for selecting the low‐rank correction matrix is very effective. Copyright © 2002 John Wiley & Sons, Ltd.     

A quadrature based method for evaluating exponential-type functions for exponential methods

We present a quadrature-based method to evaluate exponential-like operators required by different kinds of exponential integrators. The method approximates these operators by means of a quadrature formula that converges like O(e ^−cK ), with K the number of quadrature nodes, and it is useful when solving parabolic equations. The approach allows also the evaluation of the associated scalar mappings. The method is based on numerical inversion of sectorial Laplace transforms. Several numerical illustrations are provided to test the algorithm, including examples with a mass matrix and the application of the method inside the MATLAB package EXP4, an adaptive solver based on an exponential Runge–Kutta method.     

A look-ahead algorithm for the solution of general Hankel systems

Summary The solution of systems of linear equations with Hankel coefficient matrices can be computed with onlyO(n ²) arithmetic operations, as compared toO(n ³) operations for the general cases. However, the classical Hankel solvers require the nonsingularity of all leading principal submatrices of the Hankel matrix. The known extensions of these algorithms to general Hankel systems can handle only exactly singular submatrices, but not ill-conditioned ones, and hence they are numerically unstable. In this paper, a stable procedure for solving general nonsingular Hankel systems is presented, using a look-ahead technique to skip over singular or ill-conditioned submatrices. The proposed approach is based on a look-ahead variant of the nonsymmetric Lanczos process that was recently developed by Freund, Gutknecht, and Nachtigal. We first derive a somewhat more general formulation of this look-ahead Lanczos algorithm in terms of formally orthogonal polynomials, which then yields the look-ahead Hankel solver as a special case. We prove some general properties of the resulting look-ahead algorithm for formally orthogonal polynomials. These results are then utilized in the implementation of the Hankel solver. We report some numerical experiments for Hankel systems with ill-conditioned submatrices.The research of the first author was supported by DARPA via Cooperative Agreement NCC 2-387 between NASA and the Universities Space Research Association (USRA).The research of the second author was supported in part by NSF grant DRC-8412314 and Cooperative Agreement NCC 2-387 between NASA and the Universities Space Research Association (USRA).     

Fast tensor method for summation of long‐range potentials on 3D lattices with defects

In this paper, we present a method for fast summation of long‐range potentials on 3D lattices with multiple defects and having non‐rectangular geometries, based on rank‐structured tensor representations. This is a significant generalization of our recent technique for the grid‐based summation of electrostatic potentials on the rectangular L × L × L lattices by using the canonical tensor decompositions and yielding the O(L) computational complexity instead of O(L³) by traditional approaches. The resulting lattice sum is calculated as a Tucker or canonical representation whose directional vectors are assembled by the 1D summation of the generating vectors for the shifted reference tensor, once precomputed on large N × N × N representation grid in a 3D bounding box. The tensor numerical treatment of defects is performed in an algebraic way by simple summation of tensors in the canonical or Tucker formats. To diminish the considerable increase in the tensor rank of the resulting potential sum, the ?‐rank reduction procedure is applied based on the generalized reduced higher‐order SVD scheme. For the reduced higher‐order SVD approximation to a sum of canonical/Tucker tensors, we prove the stable error bounds in the relative norm in terms of discarded singular values of the side matrices. The required storage scales linearly in the 1D grid‐size, O(N), while the numerical cost is estimated by O(NL). The approach applies to a general class of kernel functions including those for the Newton, Slater, Yukawa, Lennard‐Jones, and dipole‐dipole interactions. Numerical tests confirm the efficiency of the presented tensor summation method; we demonstrate that a sum of millions of Newton kernels on a 3D lattice with defects/impurities can be computed in seconds in Matlab implementation. The tensor approach is advantageous in further functional calculus with the lattice potential sums represented on a 3D grid, like integration or differentiation, using tensor arithmetics of 1D complexity. Copyright © 2015 John Wiley & Sons, Ltd.     

Numerical computation of an analytic singular value decomposition of a matrix valued function

Summary This paper extends the singular value decomposition to a path of matricesE(t). An analytic singular value decomposition of a path of matricesE(t) is an analytic path of factorizationsE(t)=X(t)S(t)Y(t) ^T whereX(t) andY(t) are orthogonal andS(t) is diagonal. To maintain differentiability the diagonal entries ofS(t) are allowed to be either positive or negative and to appear in any order. This paper investigates existence and uniqueness of analytic SVD's and develops an algorithm for computing them. We show that a real analytic pathE(t) always admits a real analytic SVD, a full-rank, smooth pathE(t) with distinct singular values admits a smooth SVD. We derive a differential equation for the left factor, develop Euler-like and extrapolated Euler-like numerical methods for approximating an analytic SVD and prove that the Euler-like method converges.Partial support received from SFB 343, Diskrete Strukturen in der Mathematik, Universität BielefeldPartial support received from FSP Mathematisierung, Universität BielefeldPartial support received from FSP Mathematisierung, Universität BielefeldPartial support received from National Science Foundation grant CCR-8820882. Some support was also received from the University of Kansas through International Travel Fund 560478 and General Research Allocations # 3758-20-0038 and #3692-20-0038.     

Kronecker product approximations for image restoration with anti‐reflective boundary conditions

The anti‐reflective boundary condition for image restoration was recently introduced as a mathematically desirable alternative to other boundary conditions presently represented in the literature. It has been shown that, given a centrally symmetric point spread function (PSF), this boundary condition gives rise to a structured blurring matrix, a submatrix of which can be diagonalized by the discrete sine transform (DST), leading to an O(n² log n) solution algorithm for an image of size n × n. In this paper, we obtain a Kronecker product approximation of the general structured blurring matrix that arises under this boundary condition, regardless of symmetry properties of the PSF. We then demonstrate the usefulness and efficiency of our approximation in an SVD‐based restoration algorithm, the computational cost of which would otherwise be prohibitive. Copyright © 2005 John Wiley & Sons, Ltd.     

A fast solver for linear systems with displacement structure

We describe a fast solver for linear systems with reconstructible Cauchy-like structure, which requires O(rn ²) floating point operations and O(rn) memory locations, where n is the size of the matrix and r its displacement rank. The solver is based on the application of the generalized Schur algorithm to a suitable augmented matrix, under some assumptions on the knots of the Cauchy-like matrix. It includes various pivoting strategies, already discussed in the literature, and a new algorithm, which only requires reconstructibility. We have developed a software package, written in Matlab and C-MEX, which provides a robust implementation of the above method. Our package also includes solvers for Toeplitz(+Hankel)-like and Vandermonde-like linear systems, as these structures can be reduced to Cauchy-like by fast and stable transforms. Numerical experiments demonstrate the effectiveness of the software.     

On the existence of mandatory representation designs with block sizes 3 and k

A mandatory representation design MR[ν,K] is a pairwise balanced design on ν points with block sizes from K in which for each k ∈ K there is a block in the design of size k. Mendelsohn and Rees [4] investigated the existence of MR[ν,K]s, where 3 ∈ K. In this report we consider additional necessary conditions, where K = {3,k}. These conditions are proved to be sufficient for 4 ≤ k ≤ 50 with one genuine exception. © 2000 John Wiley & Sons, Inc. J Combin Designs 8: 122–131, 2000     

A fast Fourier transform based direct solver for the Helmholtz problem

This article is devoted to the efficient numerical solution of the Helmholtz equation in a two‐ or three‐dimensional (2D or 3D) rectangular domain with an absorbing boundary condition (ABC). The Helmholtz problem is discretized by standard bilinear and trilinear finite elements on an orthogonal mesh yielding a separable system of linear equations. The main key to high performance is to employ the fast Fourier transform (FFT) within a fast direct solver to solve the large separable systems. The computational complexity of the proposed FFT‐based direct solver is $\mathrm{??}(N\log N)$ operations. Numerical results for both 2D and 3D problems are presented confirming the efficiency of the method discussed.     

A direct method to solve block banded block Toeplitz systems with non-banded Toeplitz blocks

A fast solution algorithm is proposed for solving block banded block Toeplitz systems with non-banded Toeplitz blocks. The algorithm constructs the circulant transformation of a given Toeplitz system and then by means of the Sherman-Morrison-Woodbury formula transforms its inverse to an inverse of the original matrix. The block circulant matrix with Toeplitz blocks is converted to a block diagonal matrix with Toeplitz blocks, and the resulting Toeplitz systems are solved by means of a fast Toeplitz solver.The computational complexity in the case one uses fast Toeplitz solvers is equal to ξ(m,n,k)=O(mn³)+O(k³n³) flops, there are m block rows and m block columns in the matrix, n is the order of blocks, 2k+1 is the bandwidth. The validity of the approach is illustrated by numerical experiments.     

A GPU-based algorithm for efficient LES of high Reynolds number flows in heterogeneous CPU/GPU supercomputers

Αn optimized MPI+OpenACC implementation model that performs efficiently in CPU/GPU systems using large-eddy simulation is presented. The code was validated for the simulation of wave boundary-layer flows against numerical and experimental data in the literature. A direct Fast-Fourier-Transform-based solver was developed for the solution of the Poisson equation for pressure taking advantage of the periodic boundary conditions. This solver was optimized for parallel execution in CPUs and outperforms by 10 times in computational time a typical iterative preconditioned conjugate gradient solver in GPUs. In terms of parallel performance, an overlapping strategy was developed to reduce the overhead of performing MPI communications using GPUs. As a result, the weak scaling of the algorithm was improved up to 30%. Finally, a large-scale simulation (Re = 2 × 10⁵) using a grid of 4 × 10⁸ cells was executed, and the performance of the code was analyzed. The simulation was launched using up to 512 nodes (512 GPUs + 6144 CPU-cores) on one of the current top 10 supercomputers of the world (Piz Daint). A comparison of the overall computational time showed that the GPU version was 4.2 times faster than the CPU one. The parallel efficiency of this strategy (47%) is competitive compared with the state-of-the-art CPU implementations, and it has the potential to take advantage of modern supercomputing capabilities.     

Analyse d'un schéma décalé pour l'approximation de systèmes hyperboliques non conservatifs

We consider a model system made of two nonlinear equations which are non conservative. A conservation law can be obtained from these equations through linear operations only, which don't modify the shock waves. A numerical scheme based on a different mesh adapted to each variable is proposed. By choosing a shifted mesh, we have un explicit Riemann solver and we can derive a finite volume scheme. We prove a priori estimates in L^∞ norm and Total Variation for the system, which lead to a strong convergence in L¹ norm towards a solution satisfying the associated conservation law.