Structured FISTA for image restoration

In this paper, we propose an efficient numerical scheme for solving some large‐scale ill‐posed linear inverse problems arising from image restoration. In order to accelerate the computation, two different hidden structures are exploited. First, the coefficient matrix is approximated as the sum of a small number of Kronecker products. This procedure not only introduces one more level of parallelism into the computation but also enables the usage of computationally intensive matrix–matrix multiplications in the subsequent optimization procedure. We then derive the corresponding Tikhonov regularized minimization model and extend the fast iterative shrinkage‐thresholding algorithm (FISTA) to solve the resulting optimization problem. Because the matrices appearing in the Kronecker product approximation are all structured matrices (Toeplitz, Hankel, etc.), we can further exploit their fast matrix–vector multiplication algorithms at each iteration. The proposed algorithm is thus called structured FISTA (sFISTA). In particular, we show that the approximation error introduced by sFISTA is well under control and sFISTA can reach the same image restoration accuracy level as FISTA. Finally, both the theoretical complexity analysis and some numerical results are provided to demonstrate the efficiency of sFISTA.     

The generalized centro‐symmetric and least squares generalized centro‐symmetric solutions of the matrix equation AYB + CYTD = E

An n×n real matrix P is said to be a symmetric orthogonal matrix if P = P^?1 = P^T. An n × n real matrix Y is called a generalized centro‐symmetric with respect to P, if Y = PYP. It is obvious that every matrix is also a generalized centro‐symmetric matrix with respect to I. In this work by extending the conjugate gradient approach, two iterative methods are proposed for solving the linear matrix equation and the minimum Frobenius norm residual problem over the generalized centro‐symmetric Y, respectively. By the first (second) algorithm for any initial generalized centro‐symmetric matrix, a generalized centro‐symmetric solution (least squares generalized centro‐symmetric solution) can be obtained within a finite number of iterations in the absence of round‐off errors, and the least Frobenius norm generalized centro‐symmetric solution (the minimal Frobenius norm least squares generalized centro‐symmetric solution) can be derived by choosing a special kind of initial generalized centro‐symmetric matrices. We also obtain the optimal approximation generalized centro‐symmetric solution to a given generalized centro‐symmetric matrix Y₀ in the solution set of the matrix equation (minimum Frobenius norm residual problem). Finally, some numerical examples are presented to support the theoretical results of this paper. Copyright © 2011 John Wiley & Sons, Ltd.     

Generalized Stokes resolvent equations in an infinite layer with mixed boundary conditions

In this paper we prove unique solvability of the generalized Stokes resolvent equations in an infinite layer Ω₀ = ℝ^{n –1} × (–1, 1), n ≥ 2, in L^q ‐Sobolev spaces, 1 < q < ∞, with slip boundary condition of on the “upper boundary” ∂Ω⁺₀ = ℝ^{n –1} × {1} and non‐slip boundary condition on the “lower boundary” ∂Ω^–₀ = ℝ^{n –1} × {–1}. The solution operator to the Stokes system will be expressed with the aid of the solution operators of the Laplace resolvent equation and a Mikhlin multiplier operator acting on the boundary. The present result is the first step to establish an L^q ‐theory for the free boundary value problem studied by Beale [9] and Sylvester [22] in L ²‐spaces. (© 2006 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

The solvability conditions for inverse eigenproblem of symmetric and anti‐persymmetric matrices and its approximation

The problem of generating a matrix A with specified eigen‐pair, where A is a symmetric and anti‐persymmetric matrix, is presented. An existence theorem is given and proved. A general expression of such a matrix is provided. We denote the set of such matrices by ??????_Eⁿ. The optimal approximation problem associated with ??????_Eⁿ is discussed, that is: to find the nearest matrix to a given matrix A* by A∈??????_Eⁿ. The existence and uniqueness of the optimal approximation problem is proved and the expression is provided for this nearest matrix. Copyright © 2002 John Wiley & Sons, Ltd.     

Carathéodory–Julia type conditions and symmetries of boundary asymptotics for analytic functions on the unit disk

It is shown that the following conditions are equivalent for the generalized Schur class functions at a boundary point t₀ ∈ ??: 1) Carathéodory–Julia type condition of order n; 2) agreeing of asymptotics of the original function from inside and of its continuation by reflection from outside of the unit disk ?? up to order 2n + 1; 3) t₀‐isometry of the coefficients ofthe boundary asymptotics; 4) a certain structured matrix ? constructed from these coefficients being Hermitian. It is also shown that for an arbitrary analytic function, properties 2), 3), 4) are still equivalent to each other and imply 1) (© 2009 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

Least‐squares solutions of matrix inverse problem for bi‐symmetric matrices with a submatrix constraint

An n × n real matrix A = (a_ij)_{n × n} is called bi‐symmetric matrix if A is both symmetric and per‐symmetric, that is, a_ij = a_ji and a_ij = a_n+1?1,n+1?i (i, j = 1, 2,..., n). This paper is mainly concerned with finding the least‐squares bi‐symmetric solutions of matrix inverse problem AX = B with a submatrix constraint, where X and B are given matrices of suitable sizes. Moreover, in the corresponding solution set, the analytical expression of the optimal approximation solution to a given matrix A^* is derived. A direct method for finding the optimal approximation solution is described in detail, and three numerical examples are provided to show the validity of our algorithm. Copyright © 2007 John Wiley & Sons, Ltd.     

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.     

Partial Fourier approximation of the Lamé equations in axisymmetric domains

In this paper, we study the partial Fourier method for treating the Lamé equations in three‐dimensional axisymmetric domains subjected to non‐axisymmetric loads. We consider the mixed boundary value problem of the linear theory of elasticity with the displacement û , the body force f̂ ϵ (L₂)³ and homogeneous Dirichlet and Neumann boundary conditions. The partial Fourier decomposition reduces, without any error, the three‐dimensional boundary value problem to an infinite sequence of two‐dimensional boundary value problems, whose solutions û _n (n = 0, 1, 2,…) are the Fourier coefficients of û . This process of dimension reduction is described, and appropriate function spaces are given to characterize the reduced problems in two dimensions. The trace properties of these spaces on the rotational axis and some properties of the Fourier coefficients û _n are proved, which are important for further numerical treatment, e.g. by the finite‐element method. Moreover, generalized completeness relations are described for the variational equation, the stresses and the strains. The properties of the resulting system of two‐dimensional problems are characterized. Particularly, a priori estimates of the Fourier coefficients û _n and of the error of the partial Fourier approximation are given. Copyright © 1999 John Wiley & Sons, Ltd.     

The retraction algorithm for factoring banded symmetric matrices

Let A be an n × n symmetric matrix of bandwidth 2m + 1. The matrix need not be positive definite. In this paper we will present an algorithm for factoring A which preserves symmetry and the band structure and limits element growth in the factorization. With this factorization one may solve a linear system with A as the coefficient matrix and determine the inertia of A, the number of positive, negative, and zero eigenvalues of A. The algorithm requires between 1/2nm² and 5/4nm² multiplications and at most (2m + 1)n locations compared to non‐symmetric Gaussian elimination which requires between nm² and 2nm² multiplications and at most (3m + 1)n locations. Our algorithm reduces A to block diagonal form with 1 × 1 and 2 × 2 blocks on the diagonal. When pivoting for stability and subsequent transformations produce non‐zero elements outside the original band, column/row transformations are used to retract the bandwidth. To decrease the operation count and the necessary storage, we use the fact that the correction outside the band is rank‐1 and invert the process, applying the transformations that would restore the bandwidth first, followed by a modified correction. This paper contains an element growth analysis and a computational comparison with LAPACKs non‐symmetric band routines and the Snap‐back code of Irony and Toledo. Copyright © 2007 John Wiley & Sons, Ltd.     

A Robust Algorithm for Blind Total Variation Restoration

Image restoration is a fundamental problem in image processing. Blind image restoration has a great value in its practical application. However, it is not an easy problem to solve due to its complexity and difficulty. In this paper, we combine our robust algorithm for known blur operator with an alternating minimization implicit iterative scheme to deal with blind deconvolution problem, recover the image and identify the point spread function（PSF）. The only assumption needed is satisfy the practical physical sense. Numerical experiments demonstrate that this minimization algorithm is efficient and robust over a wide range of PSF and have almost the same results compared with known PSF algorithm.     

The inverse problem of bisymmetric matrices with a submatrix constraint

An n×n real matrix A is called a bisymmetric matrix if A=A^T and A=S_nAS_n, where S_n is an n×n reverse unit matrix. This paper is mainly concerned with solving the following two problems: Problem I Given n×m real matrices X and B, and an r×r real symmetric matrix A₀, find an n×n bisymmetric matrix A such that where A([1: r]) is a r×r leading principal submatrix of the matrix A. Problem II Given an n×n real matrix A^*, find an n×n matrix  in S_E such that where ∥·∥ is Frobenius norm, and S_E is the solution set of Problem I. The necessary and sufficient conditions for the existence of and the expressions for the general solutions of Problem I are given. The explicit solution, a numerical algorithm and a numerical example to Problem II are provided. Copyright © 2003 John Wiley & Sons, Ltd.     

Single‐cell compact finite‐difference discretization of order two and four for multidimensional triharmonic problems

In this article, we discuss finite‐difference methods of order two and four for the solution of two‐and three‐dimensional triharmonic equations, where the values of u,(?²u/?n²) and (?⁴u/?n⁴) are prescribed on the boundary. For 2D case, we use 9‐ and for 3D case, we use 19‐ uniform grid points and a single computational cell. We introduce new ideas to handle the boundary conditions and do not require to discretize the boundary conditions at the boundary. The Laplacian and the biharmonic of the solution are obtained as byproduct of the methods. The resulting matrix system is solved by using the appropriate block iterative methods. Computational results are provided to demonstrate the fourth‐order accuracy of the proposed methods. © 2009 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2010     

An approximation algorithm for counting contingency tables

We present a randomized approximation algorithm for counting contingency tables, m × n non‐negative integer matrices with given row sums R = (r₁,…,r_m) and column sums C = (c₁,…,c_n). We define smooth margins (R,C) in terms of the typical table and prove that for such margins the algorithm has quasi‐polynomial N^{O(ln N)} complexity, where N = r₁ + … + r_m = c₁ + … + c_n. Various classes of margins are smooth, e.g., when m = O(n), n = O(m) and the ratios between the largest and the smallest row sums as well as between the largest and the smallest column sums are strictly smaller than the golden ratio (1 + )/2 ≈? 1.618. The algorithm builds on Monte Carlo integration and sampling algorithms for log‐concave densities, the matrix scaling algorithm, the permanent approximation algorithm, and an integral representation for the number of contingency tables. © 2010 Wiley Periodicals, Inc. Random Struct. Alg., 2010     

Gauss‐Green theorem for weakly differentiable vector fields,sets of finite perimeter,and balance laws

We analyze a class of weakly differentiable vector fields F : ?ⁿ → ?ⁿ with the property that F ∈ L^∞ and div F is a (signed) Radon measure. These fields are called bounded divergence‐measure fields. The primary focus of our investigation is to introduce a suitable notion of the normal trace of any divergence‐measure field F over the boundary of an arbitrary set of finite perimeter that ensures the validity of the Gauss‐Green theorem. To achieve this, we first establish a fundamental approximation theorem which states that, given a (signed) Radon measure μ that is absolutely continuous with respect to ??^{N ? 1} on ?^N, any set of finite perimeter can be approximated by a family of sets with smooth boundary essentially from the measure‐theoretic interior of the set with respect to the measure ||μ||, the total variation measure. We employ this approximation theorem to derive the normal trace of F on the boundary of any set of finite perimeter E as the limit of the normal traces of F on the boundaries of the approximate sets with smooth boundary so that the Gauss‐Green theorem for F holds on E. With these results, we analyze the Cauchy flux that is bounded by a nonnegative Radon measure over any oriented surface (i.e., an (N ? 1)‐dimensional surface that is a part of the boundary of a set of finite perimeter) and thereby develop a general mathematical formulation of the physical principle of the balance law through the Cauchy flux. Finally, we apply this framework to the derivation of systems of balance laws with measure‐valued source terms from the formulation of the balance law. This framework also allows the recovery of Cauchy entropy flux through the Lax entropy inequality for entropy solutions of hyperbolic conservation laws. © 2008 Wiley Periodicals, Inc.     

On a quadrature algorithm for the piecewise linear wavelet collocation applied to boundary integral equations

In this paper, we consider a piecewise linear collocation method for the solution of a pseudo‐differential equation of order r=0, ?1 over a closed and smooth boundary manifold. The trial space is the space of all continuous and piecewise linear functions defined over a uniform triangular grid and the collocation points are the grid points. For the wavelet basis in the trial space we choose the three‐point hierarchical basis together with a slight modification near the boundary points of the global patches of parametrization. We choose linear combinations of Dirac delta functionals as wavelet basis in the space of test functionals. For the corresponding wavelet algorithm, we show that the parametrization can be approximated by low‐order piecewise polynomial interpolation and that the integrals in the stiffness matrix can be computed by quadrature, where the quadrature rules are composite rules of simple low‐order quadratures. The whole algorithm for the assembling of the matrix requires no more than O(N [logN]³) arithmetic operations, and the error of the collocation approximation, including the compression, the approximative parametrization, and the quadratures, is less than O(N^?(2?r)/2). Note that, in contrast to well‐known algorithms by Petersdorff, Schwab, and Schneider, only a finite degree of smoothness is required. In contrast to an algorithm of Ehrich and Rathsfeld, no multiplicative splitting of the kernel function is required. Beside the usual mapping properties of the integral operator in low order Sobolev spaces, estimates of Calderón–Zygmund type are the only assumptions on the kernel function. Copyright © 2003 John Wiley & Sons, Ltd.     

Well‐posedness of the upper convected Maxwell fluid in the limit of infinite Weissenberg number

An iteration scheme is used to show the well‐posedness of the initial‐boundary value problem for incompressible hypoelastic materials, which arise as a high Weissenberg number limit of viscoelastic fluids. We first assume that the stress is a rank‐one matrix T=qq^T, q∈?ⁿ, and develop energy estimates to show that the problem is locally well‐posed. This problem is related to incompressible ideal magnetohydrodynamics (MHD). We show that the general case T=CC^T, C∈?^n×n can be handled by a generalization of the method we developed. Copyright © 2010 John Wiley & Sons, Ltd.     

Finite element approximations to one‐phase nonlinear free boundary problem in groundwater contamination flow

In this article, we consider a single‐phase coupled nonlinear Stefan problem of the water‐head and concentration equations with nonlinear source and permeance terms and a Dirichlet boundary condition depending on the free‐boundary function. The problem is very important in subsurface contaminant transport and remediation, seawater intrusion and control, and many other applications. While a Landau type transformation is introduced to immobilize the free boundary, a transformation for the water‐head and concentration functions is defined to deal with the nonhomogeneous Dirichlet boundary condition, which depends on the free boundary function. An H¹‐finite element method for the problem is then proposed and analyzed. The existence of the approximation solution is established, and error estimates are obtained for both the semi‐discrete schemes and the fully discrete schemes. © 2006 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2006     

An implicit QR algorithm for symmetric semiseparable matrices

The QR algorithm is one of the classical methods to compute the eigendecomposition of a matrix. If it is applied on a dense n × n matrix, this algorithm requires O(n³) operations per iteration step. To reduce this complexity for a symmetric matrix to O(n), the original matrix is first reduced to tridiagonal form using orthogonal similarity transformations. In the report (Report TW360, May 2003) a reduction from a symmetric matrix into a similar semiseparable one is described. In this paper a QR algorithm to compute the eigenvalues of semiseparable matrices is designed where each iteration step requires O(n) operations. Hence, combined with the reduction to semiseparable form, the eigenvalues of symmetric matrices can be computed via intermediate semiseparable matrices, instead of tridiagonal ones. The eigenvectors of the intermediate semiseparable matrix will be computed by applying inverse iteration to this matrix. This will be achieved by using an O(n) system solver, for semiseparable matrices. A combination of the previous steps leads to an algorithm for computing the eigenvalue decompositions of semiseparable matrices. Combined with the reduction of a symmetric matrix towards semiseparable form, this algorithm can also be used to calculate the eigenvalue decomposition of symmetric matrices. The presented algorithm has the same order of complexity as the tridiagonal approach, but has larger lower order terms. Numerical experiments illustrate the complexity and the numerical accuracy of the proposed method. Copyright © 2005 John Wiley & Sons, Ltd.     

New fast divide‐and‐conquer algorithms for the symmetric tridiagonal eigenvalue problem

In this paper, two accelerated divide‐and‐conquer (ADC) algorithms are proposed for the symmetric tridiagonal eigenvalue problem, which cost O(N²r) flops in the worst case, where N is the dimension of the matrix and r is a modest number depending on the distribution of eigenvalues. Both of these algorithms use hierarchically semiseparable (HSS) matrices to approximate some intermediate eigenvector matrices, which are Cauchy‐like matrices and are off‐diagonally low‐rank. The difference of these two versions lies in using different HSS construction algorithms, one (denoted by ADC1) uses a structured low‐rank approximation method and the other (ADC2) uses a randomized HSS construction algorithm. For the ADC2 algorithm, a method is proposed to estimate the off‐diagonal rank. Numerous experiments have been carried out to show their stability and efficiency. These algorithms are implemented in parallel in a shared memory environment, and some parallel implementation details are included. Comparing the ADCs with highly optimized multithreaded libraries such as Intel MKL, we find that ADCs could be more than six times faster for some large matrices with few deflations. Copyright © 2016 John Wiley & Sons, Ltd.     

Direct algorithm for the solution of two‐sided obstacle problems with M‐matrix

A direct algorithm for the solution to the affine two‐sided obstacle problem with an M‐matrix is presented. The algorithm has the polynomial bounded computational complexity O(n³) and is more efficient than those in (Numer. Linear Algebra Appl. 2006; 13 :543–551). Copyright © 2010 John Wiley & Sons, Ltd.