Bilevel Optimization recently became a very active research area. This is mainly due to its important applications from machine learning. In this paper, we give a gentle introduction to algorithms, theory, and applications of bilevel optimization. In particular, we will discuss the history of bilevel optimization, its applications in power grid, hyper-parameter optimization, meta learning, as well as algorithms for solving bilevel optimization and their convergence properties. We will mainly discuss algorithms for solving two types of bilevel optimization problems: lower-level problem is strongly convex and lower-level problem is convex. We will discuss gradient methods and value-function-based methods. Decentralized and federated bilevel optimization will also be discussed.

Convolution-type nonlocal potentials are quite common and important in many science and engineering fields. Efficient and accurate evaluations of such potentials often bottleneck the real time simulations. The convolution kernel is usually singular or discontinuous at the origin and/or at the far field. The density is a smooth fast-decaying function, and is naturally well approximated by Fourier spectral method on a bounded rectangular domain, which is uniformly discretized in each spatial direction, with a nearly optimal complexity of O(N log N) that is inherited from the discrete Fast Fourier Transform. In some cases, there exists a strong spatial anisotropy in the density. Nonlocality, singularity and anisotropy are three challenges in convolution evaluation. The numerical problem is to compute the convolutions accurately and efficiently on such uniform mesh grid. In this article, we mainly review the state-of-art fast integral algorithms, including the NonUniform-FFT based method (NUFFT), Gaussian-Sum based method and Kernel Truncation method. All these methods achieve spectral accuracy with a FFT-like complexity O(N log N), and can be rewritten as a discrete convolution structure. The discrete convolution structure helps deal with strong anisotropy perfectly using a pair of FFT and inverse FFT(iFFT) on a twofold zero-padded density. Rigorous error estimates and extensive numerical results are shown to confirm the accuracy, efficiency and anisotropy performance.

Many practical problems in interdisciplinary sciences can be translated to the multivariable minimization problems of an energy function/functional in mathematics. There are two long-standing, critical problems in computational mathematics: finding the global minimum and finding the relationship between different minima. This paper mainly introduces the recently developed "solution landscape" concept and method. We will review the concept of solution landscape, saddle dynamics method for construction of solution landscape, and its applications on liquid crystals and quasicrystals.

In electronic structure calculations, Kohn-Sham equations rank among the most widely adopted mathematical models. However, due to the deficiency of available approximations for exchange-correlation energy, Kohn-Sham equations cannot well describe strongly correlated electrons systems at present. In recent decades, some researchers have studied the optimization models of strictly-correlated-electrons energy, starting from the strong-interaction limit of density functional theory. These models are hopeful to amend the very deficiency of Kohn-Sham equations and therefore broaden the applicability of the density functional theory. Since the curse of dimensionality resides in these models, some low-dimensional reformulations have been proposed. In this paper, we introduce the optimization models of strictly-correlated-electrons energy, highlight the research focus, and describe some lowdimensional reformulations. We also present the numerical approaches for these reformulations and shed light on the directions of future research.

The Boltzmann equation is one of the fundamental equations in kinetic theory, and serves as a basic building block connecting microscopic Newtonian mechanics and macroscopic continuum mechanics. Numerical approximation of the Boltzmann equation is a challenging problem mainly due to its high-dimensional, nonlocal, and nonlinear collision integral. Over the past 20 years, the spectral method based on Fourier series (or trigonometric polynomials) has become a popular and efficient deterministic method for solving the Boltzmann equation, manifested by its high accuracy and possibility of being accelerated by the fast Fourier transform. This paper aims to review the Fourier-Galerkin spectral method for the Boltzmann equation, stability and convergence of the method, fast algorithms, and generalizations to various Boltzmann-type collisional kinetic equations.

Many physical phenomena can be mathematically described by curvature-driven free interface motions, such as the evolution of films and foams, crystal growth, and so on. The motion of these films and interfaces often depends on their surface curvature and therefore can be described by the corresponding curvature flows and geometric evolution equations. The numerical computation and error analysis of the related free interface problems are still challenging problems in the field of computational mathematics. The parametric finite element method is a class of effective computational methods for approximating curvature flows, and it has been successful in simulating the evolution of some basic curvature flows, including mean curvature flow, Willmore flow, surface diffusion, and so on. This article focuses on the parametric finite element approximation of curvature flows-its origin, development and some current challenges.

Phase retrieval is raised in many areas, such as imaging, optics and quantum tomography etc, which attracts many attentions of experts from different areas, such as computational mathematics and data sciences etc. The aim of this paper is to introduce the basic theoretical problems in phase retrieval and also introduce many algorithms for solving phase retrieval.

Modeling the interatomic potential is one of the crucial problems in the field of molecular simulation. For a long time, the community faces the dilemma that the first-principles calculations are accurate but slow, while the empirical force fields are efficient but inaccurate. Machine learning is a promising approach to solve the dilemma because it achieves comparable accuracy with the first-principles calculations at a much lower expense. In this review, we present a general framework for developing the machine learning interatomic potentials, provide an incomplete list of recent work in this direction, and investigate the advantages and disadvantages of the reviewed approaches.

Computational neuroscience is an emerging interdiscipline and first appeared as a specific research field in the late 1980s. It is aimed to solve important scientific issues in neuroscience through mathematical modeling, theoretical analysis, and numerical simulation. On the one hand, neuroscience experiments provide the basis for the development of new mathematical models, theories, and algorithms. On the other hand, it is helpful to reveal mechanisms underlying experimental phenomena and discover new scientific laws through mathematical and quantitative analysis. The US Brain Initiative and the European Human Brain Project were both launched in 2013, while the Japan brain/minds project was launched in 2014. Recently, the China Brain Project ("One body, two wings") has also been approved by the State Council as one of the Innovation 2030 Major Science and Technology Projects. The investigations of brain and its related brain-inspired artificial intelligence are significant frontier sciences and have the leading strategic position of national competition in research and development. Because of this, computational neuroscience is regarded as a bridge between brain science and artificial intelligence and plays more and more important roles in frontier sciences and national strategic needs. The development of computational neuroscience may advance neuroscience, mathematics, physics, statistics, computer science, artificial intelligence, and other natural sciences and engineering disciplines. In addition, it can integrate the advantages of different disciplines to complement each other, and achieve important scientific breakthroughs.

In the Big Data era, with the advent of convenient automated data collection technologies, large-scale composite convex optimization problems are ubiquitous in many applications, such as massive data analysis, machine and statistical learning, image and signal processing. In this paper, we review a class of efficient proximal point algorithms for solving the large-scale composite convex optimization problems. Under the easy-to-implement stopping criteria and mild calmness conditions, we show the proximal point algorithm enjoys global and local asymptotic superlinear convergence. Meanwhile, based on the duality theory, we propose an efficient semismooth Newton method for handling the subproblems in the proximal point algorithm. Lastly, to further accelerate the proximal point algorithm, we fully exploit the nonsmooth second order information induced by the nonsmooth regularizer in the problem to achieve a dramatic reduction of the computational costs of solving the involved semismooth Newton linear systems.

The first principles electronic structure calculations have become important tools for studying the material mechanism, understanding and predicting the material properties, and have achieved great success. However, it is still full of challenge for how to design highly efficient and highly accurate computational methods to deal with larger system, how to understand the reliability and efficiency of calculation from a mathematical point of view. Based on the Kohn-Sham DFT, the key mathematical modes for electronic structure calculations are the Kohn-Sham equation or the Kohn-Sham energy functional minimization problem. In the past decades, the highly efficient algorithms design and numerical analysis have attracted the attention of many distinguished mathematicians. Our group have also focused on this field and have done several works. In this paper, we introduce recent progresses in this field, mainly about those done by our group.

This article considers the application of the weak Galerkin finite element (WG) method to linear elasticity problems. The WG method is a generalization of the traditional finite element method, which is used to solve numerical solutions of partial differential equations. In WG, the weak function, a piecewise polynomial function that is defined both inside the element and on the boundary of the element, is used as an approximate function and weak differential operators are given correspondingly. Moreover, stabilizers are introduced to keep the weak continuity of the approximate function. In the WG method, partitions could be arbitrary polygons or polyhedrons that satisfies the shape regular conditions. In addition the numerical format and the approximate function are easy to construct. In this paper, we introduce the application of the WG method in solving linear elasticity problems by solving three common problems in the numerical methods for linear elasticity problems, namely: the coerciveness, locking property, and the symmetry of stress tensor.

We live in the digital age, and data has become an essential part of our lives. Images are undoubtedly one of the most important types of data. Image inverse problems, including image denoising, deblurring, restoration, biomedical imaging, etc., are important areas in imaging science. The rapid development of computer technology has enabled us to use sophisticated mathematics and machine learning tools to design effective algorithms for image inverse problems. This paper mainly reviews three types of methods in image inverse problem, namely, applied and computational harmonic analysis method (represented by wavelets and wavelet frames), partial differential equation (PDE) method and deep learning method. We will review the modeling philosophies of these methods, explore the connections and differences among them, their advantages and disadvantages, and further discuss the feasibility and benefit of the integration of these methods.

Many problems in multi-agent systems, due to communication restrictions, need to be solved in a decentralized manner. There is no data fusion center, so we must rely on shortdistance communication between adjacent nodes to achieve the goal of the whole network. Compared with traditional (centralized) computing, decentralized computing is more suitable for distributed data, less subject to communication and computing bottlenecks, and easier to realize in some applications.
This article overviews the formulations and methods of decentralized consensus optimization. The objective of consensus optimization is that all the variables of the nodes converge to the same vector that minimizes the sum of their objective functions. This problem is solved by calculations at each node and data exchanges between adjacent nodes. Naive decentralized algorithms are much slower than their centralized counterparts. In order to make up for this gap, we review some recent methods through a unified framework of operator splitting.

The main goal of this paper is to introduce some recent developments of numerical methods for solving Hartree-Fock-like equations in the context of large basis sets. HartreeFock-like equations are an important type of equations in electronic structure theory. They appear in the Hartree-Fock theory, as well as the Kohn-Sham density functional theory with hybrid exchange-correlation functionals, and are widely used in electronic structure calculations of complex chemical and materials systems. Because of its high computational cost, Hartree-Fock-like equations are typically only used in systems consisting of tens to hundreds of electrons. From a mathematical perspective, Hartree-Fock-like equations are a system of nonlinear integro-differential equations. The computational cost is mainly due to the integral operator part, namely the Fock exchange operator. Through the development of the adaptive compression method (ACE), the projected commutator-direct inversion in the iterative subspace (PC-DⅡS) method, and the interpolative separable density fitting (ISDF) method, we demonstrate that the cost of Kohn-Sham density functional theory calculations with hybrid functionals can be significantly reduced. Using a silicon system with 1000 atoms for example, we have reduced the cost of hybrid functional calculations with a planewave basis set to a level that is close to the cost of semi-local functional calculations, which do not involve the Fock exchange operator. Meanwhile, we find that the structure of HartreeFock-like equations provides new insights for the iterative solution of one type of eigenvalue problems.

Abstract Trust region algorithms for nonlinear optimization and their convergence propertiesare discussed. Convergence results and techniques for convergence analysis are studied.An A (δ,η) descent trial step is defined, and is used to obtain a unified proof for globalconvergence of trust region algorithms

In this paper, some simple determinate conditions for nonsingular H-matrices areobtained, and some results include and generalize the related results in [1, 2, 3].

New error estimates for both moments and rotations are given for the usual Morley plateelement, which are obtained directly using nonconforming finite element techniques that differfrom Arnold-Brezzi approach. The error bound for rotations is better than that of Arnold-Brez-zi, while the bounds for moments in the two different approaches are identical.

It is shown that a second order quasi-conforming element, introduced recentlyby Han Hou-de, is nothing but the usual de Veubeke nonconforming element.

The following two kinds of inverse eigenvalue problems arising from structuraldesign are discussed. Problem SIELS: Given a real n×k matrix X_1 (1≤k

Some general conditions for nonsingularity of the interpolating matrix of the cubic sp-line are presented by using eigenvalue analysis techniques. The results have generalized andimproved those in [1], [2], [3].

First, an exact definition of sparsity and local nonlinearity for large nonlinear equations isgiven. Let ?~((1)) = {λ_(ij)~((1))} and ?~((2)) = {λ_(ij)~((2))}, where λ_(ij)~((1)) = 1 if x_j appears in f_i(x), 0 otherwise, λ_(ij)~((1)) = 1 if x_j appears as nonlinear term in f_i(x), 0 otherwise.Furthermore, E~((2)) = {(i, j)|λ_(ij)~((2)) = 1 } is defined. The quasi-newton scheme is well known x~((l+1)) = x~((l)) + αp, K~((l))p = -- F(x~((1))), l = 0, 1, 2,…,where K~((l+1)) = K~((l)) + M~((l)), M~((l))αp = r, r = F(x~((l+1)) M (α -- 1)F(x~((l))).We suppose that M~((l)) has the following form Then α_(ii)αp_i + α_(ij)αp_j = ξ_(ij)r_i, α_(ij)αp_i + α_(ij)αp_j = ξ_(ij)r_j.Minimizing under constraints ∑ j ξ_(ij) = 1, we can obtain {ξ_(ij)}, and theh A_(ij). The practical example shows that the above method has good efficiency, especially, for∑λ_(ij)~((2))<<∑λ_(ij)~((1)) << n~2.

The purpose of the present paper is to find affine intrinsic invariants of aplane Bezier curve B_n of nth order, and then to give certain geometrical interpretations forthe necessary and sufficient Conditions in order that a plane B_3 should be convex. Further,the inflexion points of a plane B_4 are disscussed with some concreteexamples.

In this paper, we generalize the results of [4] to three-dimensional vorticityequations. Two classes of schemes are constructed based on conservation and transport property.The index of generalized stability of schemes for periodic problem is estimated, which showsthe relationship of stability and conservation, solves the convergence of degenerated differenceequation. Two-level schemes are given. Computational errors of such schemes are proved. Theinterference of boundary shape, boundary condition and errors of boundary values are discussedin detail for initial-boundary value problem. Finally the dynamic relaxation method of steadyflow is proposed and its convergence is proved under certain conditions.

In order to approximate the perturbation bounds for invariant subspaces andgeneralized invariant subspaces, we suggest an iterative technique for the solution of the qua-dratic matrix equations which possesses quadratic convergence. Consequendy we obtain the cor-responding perturbation theorems. At the end we give a numerical example.

A spline in local coordinates which is more general than that in [1] isconsidered. The spline curve in such a case is independent of the reference local coordinatesystem. Hence the splines in both local coordinate in [1] and ordinary one are joined togetherin a more general sense.

This paper discusses some subspaces associated with eigenvalue and generalized eigenvalue problems. In this paper, we define“pair of generalized eigenvalue matrices”and“gene-ralized eigenmatrix”of a pair of matrices, and from these definitions we establish the concept of generalized invariant subspaces; establish the necessary and sufficient conditions for the existence and uniqueness of corresponding subspaces; give the relationship of generalized invariant subspaces with the“deflating pair”which was defined by G.W. Stewart.

We consider the finite element methods with trial functions not fulfilling essential boundary conditions. the variational formulation of“Saddle-point”type is applied. The convergent error estimates of optimal order of ccuracy are derived. The resulting schemes have same numerical stability as R tz finite element methods.

Boundary-value problems of elliptic equations may have many different mathematical formulations, equivalent in principle but not equally efficient in practice. For example, Neumann problem of Laplace equations (1),(2) is equivalent to the variational problem (3),(4). The judicious use of the latter formulation leads to the success of the FEM. The problem can also be formulated in terms of integral equations, even in many ways. They have generally the advantage of the reduction both of dimension by 1 and of the infinite domain to the finite, at the expense of increased analytical difficulty. The most wall-known reduction is the Fredholm integral equation of the second kind(5), for which w is to be solved and gives the original solution through the integral formula(6). The corresponding integral operator maps H~s(Ω)→H~s(Ω) and is, in general, not self-adjoint, so one of the characteristic and useful properties of the original problem is lost. A less-known reduction to integral equation is in the form(7), for which the boundary value u_0 of the solution u to the orginal problem is to be solved and gives u through the integral formula(8). The kernel K has the advantage of being self-adjoint and is derived from the Green’s function by double differentiation so is highly singular. It is of the type of the finite part of the divergent integral in the sense of Hadamard and maps H~s onto H~(s-1)and is thus desmoothing by orde r 1. This is advantageous rather than defective to the solution stability. Furthermore, the variational formulation equivalent to(7) is(11), (12)which can be obtained from (3),(4) through elimination of interior values of u by means of Green’s function. This form of reduction to integral equation is related to the original problem in a more natural and direct way, so it will be regarded as canonical and is more desirable in numerical approach. In fact, the idea of canonical reduction is implicitly used in FEM practice as technique of substructures. The elimination of the internal degrees of freedom is precisely a discrete analog of the canonical reduction and the resulted algebraic system containing solely the boundary degrees of freedom is precisely a discrete analog of the Hadamard integral equation. Recently, an elegant scheme of infinite similar elements has been proposed for the solution of crack singularity problems. They are equally well suited for concave corners, intersection of several interfaces, infinite domains and also the usual closed domain of regularity. For all these cases, it can be shown that, under certain uniformity condition, a conforming finite element in infinite similar triangulation converges with its nominal order of accuracy without deterioration. This elimination of infinite number of the interior degrees of freedom is another example of discrete analog of the canonical reduction. Fig. 1 affords an example problem containing various kinds of singularity and infinite domain. It can be grossly divided into 5 substructures using infinite triangulation for each. This suggests an economy of problem preparation, storage space and volume of computation. Fig. 2 is an infinite triangulation of the unit circle, the finite algebraic system for the boundary unknowns after the elimination of infinite many interior unknowns gives a discrete analog of the Hadamard integral equation(14) with the finite part kernel 1/sin~2θfor the unit circle. This is an example of solving integral equation without explicit use of integral equation, also that of treating finite parts without explicit presentation of finite parts.

In this note, the sufficient conditions for the existence of bifurcation solutions for boundary valued problem u(-R)=u(R)=u_0 of ordinary differential equations (K(u)u′)′+λF(u)=0 are considered. Several numerical examples are given.

The theory of finite elements has been established since the early sixties and hasbeen developed to a certain degree of completeness and sophistication for the classicalcontinuous (conforming) cose. The theory for the discontinuous (nonconforming)case is still in a less satisfactory state, although important progress has been made.The present work deals with the theoretical foundation of the discontinuous finiteelements. In section 1, Poincare inequalities for discontinuous functions are given. Theydiffer from the classical ones by an additional term of the integral of jump valuessquared with a constant which measures the density of distribution of discontinuities.On this basis, in sections 2 and 3, injection theorems--discrete analogs of the classicalones--for the discontinuous finite element functions spaces can be established for thecase of formal Sobolev norm (discontinuity discarded) as well as for the case of normcentaining additional penalty (counting the discontinuity). For the first case, a cer-tain condition of weak discontinuity is imposed and this condition is satisfied practi-cally by all the non-conforming elements now in use. In the second case, the condi-tion of weak discontinuity may be violated, i.e., the discontinuity may be arbitrarilystrong. This suggests, among others. two kinds of policy for using discontinuouselements: the policy of tolerance--this is the usual method--in ease of weak dis-continuity and the policy of suppression--this is the penalty method--in case ofstrong discontinuity. In section 4 a general convergence theorem of the penalty method for solvingelliptic equations of order 2m is given to the effect that it is always convergentwhen the finite element interpolation operator is exact to the degree k≥m and thepenalty parameters p_i satisfy 1<2k-2i+1, i = 0, 1,..., m-1. The choicep_i= k+1-i gives the best order O (h(k-m+1)/2) of convergence for sufficiently smoo-th solutions. As a result, all membrane elements of degree of accuracy k≥2 canbe used as plate elements with convergence order O(h(k-1)/2).

In the usual theory of partial differential equations,one considers in R~n a domain of homogeneous dimension n,on which the differential equations are defined,and the boundary of dimension n-1,on which the boundary conditions are prescribed.It is natural and desirable to extend such a setting to the case where the domain is of heterogeneous dimensions,i.e.,it consists of a finite number of pieces of different dimensions,suitably connected together,with differential equations on each piece cou- pled through incidence relations and eventually supplemented by boundary conditions on the remaining boundaries.This is actually the mathematical situation in many engineering problems,and the great geometrical and analytical complexities herein encountered should be trickled in a proper mathematical way. In section 1 we define a composite manifold as a closed complex of cells of dif- ferent dimensions,each cell being a connected smooth orientable Riemannian mani- fold with piecewisc smooth boundary,i.e.,the boundary consists of a finite number of cells of 1 dimension lower.In a composite manifold,a subcomplex is defined as a composite structure Ω when its closure coincides with the underlying composite manifold and certain strong connectedness property is satisfied;and another subcom- plex is suitably defined as the boundary Ω of the composite structure Ω.The cou- pled differential equations will be defined on each cell of Ω and the boundary condi- tions will be prescribed on each cell of Ω. In section,2,a product of Sobolev spaces of order 1 corresponding to all the cells of a composite structure is introduced and a closed subspacc is specified by cer- tain link condition.For this subspace,some injection theorems in sense of Sobolev can be established and a standard elliptic variational problem is introduced and leads to a system of coupled Poissou equations on a conlpositc manifold in R~n which is a natural extension of the classical Poisson equation and is applicable to the heat transfer,diffusion on complex structures. In section 3,considerations analogous to those of section 2 lead to another product of Sobolev spaces for a composite structure in R~3 and the corresponding injection theorems.This gives a precise mathematical foundation for the composite elastic structures. Differential equations on composite manifolds seem to have wide applications. Some relevent theoretial problems worthy of further study are indicated.

Since the assumption of strang's uniformity on the basis functions is replaced by the uniformity assumption on interpolated nodes, the opthnal approximablity error estimates are derived. Considering an extreme value of certain set function, the "inverse relations" of pieeewise polynomials also are achieved generally.

A mathematical theory of thick plate is established on another mechanical hypothesis differed from a thin plate. In this paper a statical analysis for thick plate is given. For an arbitrary plate of polygonal shape with simply supported edges it is shown that the solution of thick plate can be "expressed explicitly by that of thin plate. It is then found that the shear forces of thick plate in the Reissner model are identical with those of thin plate, while in the washizu model of thick plate the shear forces as well as bending moments and angular rotations all are same as the classical thin plate tbeory.

In this paper a finite element method based on cubic B-spline is presented to obtainnig approximate solution for the equilibrium problems of elastic composite structures on regular regions.A computational scheme well suited for various types of boundary conditions is derived.It is easily carried out by a computer.In comparison with the usual nodal finite element method, the main features of the present method are higher accuracy and more economy in computing storage and time requirements.The effect of imposing natural boundary condition as constraints is considered and it is shown by a plate bending problem that this effect is quite superficial.Several numerical results are presented.

This article is intended to give a new iterative method to calculate the combined stiffness matrix in the infinite element method.The convergence of this iterative method is proved.Three simple examples are given at the end of this article.They are quite fast in convergence rate and fairly high in precision.