A Preconditional Manifold Accelerated Optimization Method for Spherically Constrained Energy Functions

Abstract

This invention discloses a preconditional manifold acceleration optimization method for spherically constrained energy functionals, belonging to the field of scientific computing and numerical optimization. The method performs numerical discretization on the target energy functional, constructing a discrete-form spherically constrained optimization problem; defines the geometric structure of the unit spherical manifold and its corresponding tangent space, and determines the orthogonal projection rules of the tangent space; introduces an auxiliary scalar sequence to control the momentum coefficient, constructs prediction points through a manifold shrinkage operator, calculates the Euclidean gradient at the prediction point and projects it to obtain the Riemann gradient, and generates a preconditional search direction after correction by a preconditional operator; performs a line search along the search direction to determine the step size, and generates the next iteration point through manifold shrinkage mapping. This invention can integrate momentum acceleration and preconditioning mechanisms while strictly maintaining spherical constraints, improving convergence speed and adaptability to ill-conditioned problems, reducing the solution overhead of large-scale discrete systems, and adapting to various minimization scenarios of energy functionals with unit modulus constraints.

Description

Technical Field

[0001] This invention relates to the fields of scientific computing and numerical optimization, and in particular to a preconditional manifold acceleration optimization method for spherically constrained energy functionals. Background Technology

[0002] The energy functional minimization problem under spherical constraints is a fundamental problem in scientific computing and engineering optimization, widely applicable in various scenarios such as quantum physics simulation, electronic structure computation, large-scale matrix eigenvalue solving, signal processing, and beam optimization. For this type of constrained optimization problem, several mature solution techniques have emerged, including gradient-based methods, conjugate gradient methods, and Newton and quasi-Newton methods based on Riemannian geometry. Related theoretical research continues to advance, and various improved schemes are constantly emerging, adaptable to solution requirements with different dimensions and degrees of nonlinearity. With the expansion of computational scale and the increasing demands for solution accuracy, the application demand for this type of optimization method in large-scale discrete systems continues to grow. Related technologies are also constantly developing towards lower computational overhead and higher parallel adaptability, providing diversified technical support for the numerical solution of various engineering and scientific problems.

[0003] Existing solution techniques still have several shortcomings in practical applications. First, some methods rely on projection or normalization operations after iteration to maintain constraint feasibility when dealing with spherical constraints, which can easily affect the continuity of the search direction and make it difficult to achieve a good balance between constraint preservation and directional effectiveness. Second, conventional first-order manifold optimization methods mostly only utilize the local gradient information of the current iteration point and fail to fully combine the displacement and gradient information of historical iterations to build an acceleration mechanism, resulting in significant room for improvement in convergence speed. Third, for ill-conditioned systems with high condition numbers after discretization, existing methods fail to deeply integrate effective preconditioning strategies with the manifold geometric framework, making it difficult to effectively improve the spectral distribution characteristics of the objective functional and prone to convergence stagnation in large-scale discrete scenarios. In addition, most existing techniques only improve a single performance dimension and lack a unified optimization framework that can simultaneously take into account manifold constraint preservation, momentum acceleration convergence, and preconditioning to improve ill-conditionedness, making it difficult to fully adapt to the solution requirements of high-dimensional non-convex constraint optimization. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the prior art and provide a preconditional manifold acceleration optimization method for spherically constrained energy functionals.

[0005] The objective of this invention is achieved through the following technical solution: A preconditional manifold acceleration optimization method for spherically constrained energy functionals is provided, which includes the following steps: S1. Perform numerical discretization on the target energy functional. The numerical discretization process includes finite difference processing, spectral method processing and finite element processing to generate a discrete energy function and construct a discrete form of spherical constrained optimization problem. S2. For the spherical constraint optimization problem, define the geometric structure of the unit spherical manifold and the corresponding tangent space of the unit spherical manifold. The tangent space is the set of all vectors orthogonal to the current iteration point. Determine the computation rules of the orthogonal projection operator of the tangent space. S3. Introduce an auxiliary scalar sequence to control the momentum coefficient, construct prediction points through manifold shrinkage operators, calculate the Euclidean gradient at the prediction points based on discrete energy functions, project the Euclidean gradient to the tangent space through the tangent space orthogonal projection operator to obtain the Riemann gradient, and modify the Riemann gradient based on the preconditioning operator to generate the preconditioning search direction. S4. Perform a line search along the precondition search direction to determine the iteration step size. Map the update vector to the unit spherical manifold using the manifold shrinking operator to generate the next iteration point. The manifold shrinking operator satisfies the first-order approximation property.

[0006] Furthermore, constructing the prediction points in step S3 includes the following sub-steps: S31. Update the value of the auxiliary scalar sequence according to the Nesterov acceleration rule. The auxiliary scalar sequence is used to control the degree of momentum decay and generate the momentum coefficient for the corresponding iteration step. S32. The momentum term is generated by linearly combining the momentum coefficient with the displacement vector generated in the previous iteration. S33. The momentum term is mapped from the tangent space to the unit spherical manifold by the manifold shrinking operator, generating the prediction point for the corresponding iteration step.

[0007] Furthermore, calculating the Riemann gradient in step S3 includes the following sub-steps: S301. Read the expression of the discrete energy function and calculate the Euclidean gradient of the discrete energy function at the prediction point. The Euclidean gradient is the derivative of the discrete energy function in Euclidean space. S302. Obtain the tangent space orthogonal projection operator corresponding to the prediction point, and project the Euclidean gradient to the tangent space corresponding to the prediction point through the tangent space orthogonal projection operator; S303. Generate the Riemann gradient at the prediction point. The Riemann gradient is the gradient of the target energy functional on the unit spherical manifold.

[0008] Furthermore, the generation of preconditional search directions in step S3 includes the following sub-steps: S311. Construct a positive definite preconditioning operator. The preconditioning operator is an approximation operator of the Riemann-Hessian matrix. The preconditioning operator satisfies the positive definiteness requirement and has a structure that supports fast inversion calculation. S312. Input the Riemann gradient into the preconditioning operator and perform a linear transformation on the Riemann gradient through the preconditioning operator; S313. Generate a modified precondition search direction, which is the descent direction in the tangent space.

[0009] Furthermore, generating the next iteration point in step S4 includes the following sub-steps: S41. Perform a line search along the precondition search direction. The line search includes Armijo line search and Wolfe line search. Determine the iteration step size that satisfies the set conditions. S42. Perform a scalar multiplication operation between the iteration step size and the precondition search direction to generate the update vector for the corresponding iteration step; S43. The update vector is mapped from the tangent space of the current iteration point to the unit spherical manifold using the manifold shrinking operator to generate the next iteration point.

[0010] Furthermore, in step S1, the spherical constraint optimization problem includes a quantum physics computation problem, an electronic structure computation problem, a numerical algebraic eigenvalue problem, and a signal processing constraint optimization problem. The target energy functional corresponding to the spherical constraint optimization problem all contain a unit modulus constraint condition, which requires that the Euclidean norm of the optimization variable equals a set value. The target energy functional corresponding to the quantum physics computation problem contains a kinetic energy term, an external potential term, and a nonlinear interaction term. The target energy functional corresponding to the electronic structure computation problem contains a single-electron kinetic energy term, an external potential term, and an electron interaction term. The target energy functional corresponding to the numerical algebraic eigenvalue problem is in the form of a matrix Rayleigh quotient. The target energy functional corresponding to the signal processing constraint optimization problem is in the form of a signal separation cost function.

[0011] Furthermore, in step S31, when the inner product of the momentum direction and the current gradient direction is greater than a set threshold, the value of the momentum term is reset to zero, and the value of the auxiliary scalar sequence is reset to the initial value; the momentum direction is the direction after projection of the search direction of the previous iteration step into the tangent space, and the current gradient direction is the Riemann gradient direction at the prediction point; the inner product value of the momentum direction and the current gradient direction is calculated, and the inner product value is compared with the set threshold; when the inner product value is greater than the set threshold, the current momentum accumulation process is terminated, and the momentum update process is restarted; when the inner product value is less than or equal to the set threshold, the current values ​​of the momentum term and the auxiliary scalar sequence are maintained, and the momentum accumulation operation continues.

[0012] Furthermore, in step S311, the preconditioning operator includes a diagonal dominant operator, an incomplete Cholesky decomposition operator, and a differential principal part approximation operator. The diagonal dominant operator extracts all elements on the diagonal of the Riemann-Hessian matrix to form a diagonal matrix, and uses the inverse of the diagonal matrix as the preconditioning operator. The incomplete Cholesky decomposition operator performs a sparse decomposition operation on the Riemann-Hessian matrix, retains non-zero elements with a set sparsity, generates a lower triangular decomposition matrix, and constructs the preconditioning path based on the lower triangular decomposition matrix. The differential principal part approximation operator extracts the second-order differential principal part operator corresponding to the energy functional and uses the Fourier transform method to realize the fast inversion calculation of the differential principal part operator.

[0013] Furthermore, in step S41, the line search uses the Barzilai-Borwein step size rule to determine the iteration step size, and calculates the iteration step size using the displacement vector and gradient difference vector of the previous iteration. The Barzilai-Borwein step size rule includes two calculation forms. The first calculation form uses the ratio of the inner product of the displacement vectors to the inner product of the displacement vector and the gradient difference vector as the iteration step size. The second calculation form uses the ratio of the inner product of the displacement vector and the gradient difference vector to the inner product of the gradient difference vector as the iteration step size. The displacement vector is calculated from the position difference between the current iteration point and the previous iteration point, and the gradient difference vector is calculated from the difference between the Riemann gradient of the current iteration point and the Riemann gradient of the previous iteration point. The line search process selects one of the calculation forms to generate the iteration step size.

[0014] Furthermore, in step S4, the manifold shrinkage operator includes a normalization shrinkage operator and a Riemann exponent mapping operator. The execution flow of the normalization shrinkage operator is as follows: the current iteration point is added to the tangent space update vector to obtain an intermediate vector, the intermediate vector is normalized using the Euclidean norm, and the normalized vector is output as the mapping result. The execution flow of the Riemann exponent mapping operator is as follows: the manifold is moved by a set length along the geodesic direction passing through the current iteration point on the unit spherical manifold, and the endpoint of the movement is output as the mapping result. Both the normalization shrinkage operator and the Riemann exponent mapping operator satisfy the first-order approximation property and can map the tangent space vector to the surface of the unit spherical manifold.

[0015] The beneficial effects of this invention are: (1) The Riemannian manifold geometry framework is used to handle the spherical constraint optimization problem. Tangent space projection and manifold shrinkage operation are used to maintain the constraint feasibility of the iteration point, eliminating the need for additional normalization correction steps and improving the numerical stability of the iteration process. (2) The preconditioning operator is applied to the Riemann gradient to complete the direction correction, optimize the spectral distribution characteristics of the objective functional, alleviate the ill-conditioned problem of large-scale discrete systems, and effectively reduce the number of iterations for high-precision solutions; (3) The Nesterov momentum mechanism is integrated to construct prediction points to accelerate iteration, while being compatible with various line search and shrinkage operator forms, and can cover multiple energy functional solution scenarios with unit modulus constraints. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the steps of a preconditional manifold acceleration optimization method for spherically constrained energy functionals. Figure 2 The flowchart of the preconditional manifold acceleration optimization method for spherically constrained energy functionals provided in the embodiment is shown. Detailed Implementation

[0017] The technical solution of the present invention will be clearly and completely described below with reference to the embodiments. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] Example 1 See Figure 1 This embodiment provides a preconditional manifold acceleration optimization method for spherically constrained energy functionals, which includes the following steps: S1. Perform numerical discretization on the target energy functional. The numerical discretization process includes finite difference processing, spectral method processing and finite element processing to generate a discrete energy function and construct a discrete form of spherical constrained optimization problem. S2. For the spherical constraint optimization problem, define the geometric structure of the unit spherical manifold and the corresponding tangent space of the unit spherical manifold. The tangent space is the set of all vectors orthogonal to the current iteration point. Determine the computation rules of the orthogonal projection operator of the tangent space. S3. Introduce an auxiliary scalar sequence to control the momentum coefficient, construct prediction points through manifold shrinkage operators, calculate the Euclidean gradient at the prediction points based on discrete energy functions, project the Euclidean gradient to the tangent space through the tangent space orthogonal projection operator to obtain the Riemann gradient, and modify the Riemann gradient based on the preconditioning operator to generate the preconditioning search direction. S4. Perform a line search along the precondition search direction to determine the iteration step size. Map the update vector to the unit spherical manifold using the manifold shrinking operator to generate the next iteration point. The manifold shrinking operator satisfies the first-order approximation property.

[0019] In some embodiments, constructing the prediction point in step S3 includes the following sub-steps: S31. Update the value of the auxiliary scalar sequence according to the Nesterov acceleration rule. The auxiliary scalar sequence is used to control the degree of momentum decay and generate the momentum coefficient for the corresponding iteration step. S32. The momentum term is generated by linearly combining the momentum coefficient with the displacement vector generated in the previous iteration. S33. The momentum term is mapped from the tangent space to the unit spherical manifold using the manifold shrinking operator, generating the prediction point for the corresponding iteration step.

[0020] In some embodiments, calculating the Riemann gradient in step S3 includes the following sub-steps: S301. Read the expression of the discrete energy function and calculate the Euclidean gradient of the discrete energy function at the prediction point. The Euclidean gradient is the derivative of the discrete energy function in Euclidean space. S302. Obtain the tangent space orthogonal projection operator corresponding to the prediction point, and project the Euclidean gradient to the tangent space corresponding to the prediction point through the tangent space orthogonal projection operator; S303. Generate the Riemann gradient at the prediction point. The Riemann gradient is the gradient of the target energy functional on the unit spherical manifold.

[0021] In some embodiments, generating the precondition search direction in step S3 includes the following sub-steps: S311. Construct a positive definite preconditioning operator. The preconditioning operator is an approximation operator of the Riemann-Hessian matrix. The preconditioning operator satisfies the positive definiteness requirement and has a structure that supports fast inversion calculation. S312. Input the Riemann gradient into the preconditioning operator and perform a linear transformation on the Riemann gradient through the preconditioning operator; S313. Generate a modified precondition search direction, which is the descent direction in the tangent space.

[0022] In some embodiments, generating the next iteration point in step S4 includes the following sub-steps: S41. Perform a line search along the precondition search direction. The line search includes Armijo line search and Wolfe line search. Determine the iteration step size that satisfies the set conditions. S42. Perform a scalar multiplication operation between the iteration step size and the precondition search direction to generate the update vector for the corresponding iteration step; S43. The update vector is mapped from the tangent space of the current iteration point to the unit spherical manifold using the manifold shrinking operator to generate the next iteration point.

[0023] In some embodiments, in step S1, the spherical constraint optimization problem includes a quantum physics computation problem, an electronic structure computation problem, a numerical algebraic eigenvalue problem, and a signal processing constraint optimization problem; the target energy functional corresponding to the spherical constraint optimization problem all contain a unit modulus constraint condition, which requires that the Euclidean norm of the optimization variable equals a set value; the target energy functional corresponding to the quantum physics computation problem contains a kinetic energy term, an external potential term, and a nonlinear interaction term; the target energy functional corresponding to the electronic structure computation problem contains a single-electron kinetic energy term, an external potential term, and an electron interaction term; the target energy functional corresponding to the numerical algebraic eigenvalue problem is in the form of a matrix Rayleigh quotient; and the target energy functional corresponding to the signal processing constraint optimization problem is in the form of a signal separation cost function.

[0024] In some embodiments, in step S31, when the inner product of the momentum direction and the current gradient direction is greater than a set threshold, the value of the momentum term is reset to zero, and the value of the auxiliary scalar sequence is reset to the initial value; the momentum direction is the direction after projection of the search direction of the previous iteration step into the tangent space, and the current gradient direction is the Riemann gradient direction at the prediction point; the value of the inner product of the momentum direction and the current gradient direction is calculated, and the value of the inner product is compared with the set threshold; when the value of the inner product is greater than the set threshold, the current momentum accumulation process is terminated, and the momentum update process is restarted; when the value of the inner product is less than or equal to the set threshold, the current values ​​of the momentum term and the auxiliary scalar sequence are maintained, and the momentum accumulation operation continues.

[0025] In some embodiments, in step S311, the preconditioning operator includes a diagonal dominant operator, an incomplete Cholesky decomposition operator, and a principal part approximation operator. The diagonal dominant operator extracts all elements on the diagonal of the Riemann-Hessian matrix to form a diagonal matrix, and uses the inverse of the diagonal matrix as the preconditioning operator. The incomplete Cholesky decomposition operator performs a sparse decomposition operation on the Riemann-Hessian matrix, retains non-zero elements with a set sparsity, generates a lower triangular decomposition matrix, and constructs the preconditioning path based on the lower triangular decomposition matrix. The principal part approximation operator extracts the second-order principal part operator corresponding to the energy functional and uses the Fourier transform method to realize the fast inversion calculation of the principal part operator.

[0026] In some embodiments, in step S41, the line search uses the Barzilai-Borwein step size rule to determine the iteration step size, and calculates the iteration step size using the displacement vector and gradient difference vector of the previous iteration. The Barzilai-Borwein step size rule includes two calculation forms. The first calculation form uses the ratio of the inner product of the displacement vectors to the inner product of the displacement vector and the gradient difference vector as the iteration step size. The second calculation form uses the ratio of the inner product of the displacement vector and the gradient difference vector to the inner product of the gradient difference vector as the iteration step size. The displacement vector is calculated from the position difference between the current iteration point and the previous iteration point, and the gradient difference vector is calculated from the difference between the Riemann gradient of the current iteration point and the Riemann gradient of the previous iteration point. The line search process selects one of the calculation forms to generate the iteration step size.

[0027] In some embodiments, in step S4, the manifold shrinkage operator includes a normalization shrinkage operator and a Riemann exponential mapping operator. The normalization shrinkage operator is executed by adding the current iteration point to the tangent space update vector to obtain an intermediate vector, performing Euclidean norm normalization on the intermediate vector, and outputting the normalized vector as the mapping result. The Riemann exponential mapping operator is executed by moving a set length along the geodesic direction passing through the current iteration point on the unit spherical manifold, and outputting the endpoint position as the mapping result. Both the normalization shrinkage operator and the Riemann exponential mapping operator satisfy the first-order approximation property and can map the tangent space vector to the surface of the unit spherical manifold.

[0028] Example 2 This embodiment provides a specific implementation process for a preconditional manifold acceleration optimization method for spherically constrained energy functionals. This implementation process is based on the Riemannian manifold geometric framework and integrates momentum acceleration mechanisms and preconditioning techniques. It can iteratively solve the minimization problem of energy functionals with unity modulus constraints. The specific implementation process is as follows: Step 1. Constructing the Discrete Constraint Problem: Step 1.1. Continuous Problems and Discrete Processing: The continuous form of the spherically constrained energy functional minimization problem aims to minimize the energy functional while requiring the optimization variables to satisfy the unit modulus constraint. In some specific implementations, the continuous problem has a unified mathematical expression, representing an abstract generalization of various problems in scientific computing and engineering. Before solving, it is necessary to clarify the specific form of the target energy functional and the constraint range. The calculation formula for the continuous form of the target problem is as follows: ; in, Let be the target energy functional, representing the continuous energy expression to be minimized; Let be the optimization variable, representing the solution function in continuous space. This problem corresponds to a minimization model with spherical constraints, whereby the optimization variables must lie on a unit spherical manifold. Discretization is a necessary prerequisite for numerical solution.

[0029] In some specific implementations, the constraint set of a continuous problem can be defined in a unified form. The constraint set limits the feasible region of the optimization variables, and all solution processes must be performed within the feasible region. The formula for calculating the constraint set is: ; in, Let be a unit spherical manifold, representing the set of all optimization variables that satisfy the unit modulus constraint; The Euclidean norm is used to measure the magnitude of a vector. The equation indicates that the square of the Euclidean norm of the optimization variable equals a set value, corresponding to the geometric definition of spherical constraints. This set of constraints forms the geometric basis of the manifold optimization framework, and all iterative operations must ensure that the iteration point falls within this set.

[0030] In some specific implementations, numerical discretization methods are used to spatially discretize the continuous energy functional. Numerical discretization processes include finite difference processing, spectral methods, and finite element methods. After processing, a discrete energy function is generated, from which a finite-dimensional discrete-constraint optimization problem is constructed. The calculation formula for the discrete-constraint optimization problem is as follows: ; in, Let be the discrete energy function, representing the finite-dimensional energy expression obtained after discretization. For discrete optimization variables, represents the solution vector corresponding to the discrete system; The degrees of freedom of the discrete system represent the total number of independent variables contained in the discrete system. Let N be an N-dimensional complex vector space, representing the numerical space in which the discrete variables reside; if the problem to be solved is in the real number domain, then the corresponding space is replaced by an N-dimensional real vector space. The discrete form of the constrained optimization problem is the solution object of the subsequent manifold optimization framework. All iterative operations are performed within the space formed by the discrete variables, and the discretization precision can be adjusted according to the actual solution requirements.

[0031] Step 2. Definition of manifold geometry: Step 2.1. Definition of a unit spherical manifold: The constraint set is defined as a unit spherical manifold, where all feasible iteration points lie on its surface. The manifold structure determines the feasible region of the optimization problem and forms the geometric basis for subsequent tangent space projection and shrinkage operations. A Riemannian manifold is a smooth manifold equipped with Riemannian metrics. In this embodiment, a unit spherical manifold is used as the constrained feasible region, transforming the constrained optimization problem into an unconstrained optimization problem on the manifold. All iterative operations follow the geometric rules of the manifold. In some specific implementations, the unit spherical manifold in discrete dimensions has a clear set expression, and any point on the manifold satisfies the unit modulus constraint. The geometric characteristics of the manifold determine the orthogonality of the tangent space. All iterative update operations must ensure that the final result falls on the manifold surface to avoid iteration points deviating from the feasible region. The formula for calculating the unit spherical manifold is: ; in, Let be a unit spherical manifold, representing the set of all optimization variables that satisfy the constraints; Let be a point on the manifold, representing a discrete optimization variable satisfying the unit modulus constraint. The unit spherical manifold transforms the constrained optimization problem into an unconstrained optimization problem on the manifold, eliminating the need to introduce additional penalty functions or Lagrange multipliers to handle constraints, thus simplifying the computational steps of the iterative process.

[0032] Step 2.2. Definition of tangent space and projection operators: Let the current iteration point be a point on the manifold. Its corresponding tangent space is defined as the set of all vectors orthogonal to the current iteration point. The tangent space is a linear approximation space of the manifold at the current iteration point, and all gradient calculations and direction updates are performed within the tangent space. The tangent space is the set of vectors corresponding to the tangent plane of the manifold at a given point. In this embodiment, it is used to hold vector variables such as gradient, momentum, and search direction, ensuring that the direction update conforms to the local geometric properties of the manifold. In some specific implementations, any vector in the tangent space is orthogonal to the normal vector at that point. Gradient projection, momentum update, and search direction construction must all be completed within the tangent space to ensure that the direction update conforms to the geometric constraints of the manifold. The formula for calculating the tangent space is: ; in, For iteration points The tangent space at a point represents the manifold. At point The set of all tangent vectors at a given point; Let be a vector in the tangent space, representing any directional vector within the tangent plane; This represents the standard complex inner product; if it is in the real number field, it corresponds to the standard dot product operation. This indicates the real part operation, used to calculate the real part of the complex inner product; Let be the iteration point of the nth iteration, representing the reference point on the manifold at the current step; the equation indicates that the real part of the inner product of the tangent vector and the iteration point is zero, corresponding to the geometric relationship of orthogonal vectors. The tangent space is the bridge connecting the gradient in Euclidean space and the gradient in the manifold. Through the projection operation, the derivative information in Euclidean space can be converted into the descent direction on the manifold.

[0033] In some specific implementations, the orthogonal projection operator is used to project any vector onto the tangent space. It is the core operation for calculating the Riemann gradient. The projection operation removes the normal component of the vector while retaining the tangent component, ensuring that the output vector lies within the tangent space. The formula for calculating the orthogonal projection operator is: ; in, For point The orthogonal projection operator at point represents projecting a vector onto the point. The operational rules corresponding to the tangent space; Let be the vector to be projected, representing any vector that needs to be projected into the tangent space; Let be the reference point on the manifold, and let represent the reference point corresponding to the tangent space. The orthogonal projection operator can be implemented through simple algebraic operations, without the need for complex matrix decomposition operations, thus reducing the computational cost of single-step iteration.

[0034] Step 3. Generate preconditional acceleration direction: Step 3.1. Prediction point construction processing: An auxiliary scalar sequence is introduced to control momentum decay. Initial values ​​for the auxiliary scalar sequence are initialized. For each iteration, momentum coefficient calculation, momentum term generation, and prediction point mapping are performed sequentially. The basic Riemann gradient method iteratively updates along the negative Riemann gradient direction, utilizing only the local gradient information of the current point, resulting in a slower convergence speed. Its iterative update calculation formula is: ; in, This is the updated iteration point; This is the shrinkage operator at the current iteration point; This is the iteration step size; Let be the Riemann gradient at the current iteration point. The Riemann gradient method is a fundamental first-order method for manifold optimization. It has a simple structure but limited convergence speed, making it difficult to adapt to the needs of large-scale, high-precision solutions.

[0035] The Riemann conjugate gradient method adds a historical direction transport term to the gradient direction in an attempt to improve convergence speed, but it is susceptible to constraint projection that violates conjugacy. Its search direction calculation formula is: ; in, This indicates the search direction for the current step. is the conjugate parameter used to control the weight of the historical direction; This is a vector transfer operator used to transfer historical search directions to the current tangent space; This serves as the search direction for the previous step. The conjugate property of the Riemann conjugate gradient method is easily destroyed by nonlinear shrinkage operations, and convergence stagnation is likely to occur in complex nonconvex problems. In this embodiment, the Nesterov momentum mechanism is used to replace the construction of the conjugate direction, which can avoid the stability problems caused by the destruction of conjugate property.

[0036] In some specific implementations, the momentum coefficient is updated according to the Nesterov acceleration rule. The iterative update of the auxiliary scalar sequence adopts a fixed recursive formula, calculating the scalar value of the current step from the scalar value of the previous step, and then deriving the momentum coefficient from the scalar values ​​of the two adjacent steps. The momentum coefficient is used to control the degree of influence of historical iteration information on the current direction. The calculation formulas for the auxiliary scalar sequence and the momentum coefficient are as follows: ; in, This is an auxiliary scalar corresponding to the nth iteration, used to recursively generate the momentum coefficient; The auxiliary scalar corresponding to the (n+1)th iteration is obtained recursively from the auxiliary scalar of the previous iteration; This is the momentum coefficient corresponding to the nth iteration, used to control the weight of the momentum term. The Nesterov acceleration rule, by introducing inertia weights to accumulate historical iteration information, can accelerate the convergence speed of the iteration process and reduce the number of iterations required to reach the convergence condition.

[0037] In some specific implementations, the momentum term is obtained by projecting the difference between two adjacent iteration points into the tangent space. The momentum term carries the displacement information from previous iterations and is the carrier for momentum acceleration. In the initial iteration stage, when there is no preceding displacement information, the momentum term is zero. The formula for calculating the momentum term is: ; in, For the momentum term corresponding to the nth iteration, it represents the tangent vector after the historical displacement information is projected onto the tangent space; Let be the iteration point of the nth iteration, representing a point on the manifold at the current step; Let be the iteration point of the (n-1)th iteration, representing a point on the manifold of the previous step; For iteration points The orthogonal projection operator at the point is used to project the displacement difference onto the tangent space. When there is no previous iteration point in the initial iteration, the previous iteration point is set to be the same as the initial iteration point. At this time, the momentum term is zero to ensure the stability of the initial stage of the iteration.

[0038] In some specific implementations, the prediction point is obtained by mapping the weighted momentum term back to the manifold using a manifold shrinking operator. The prediction point incorporates inertial information from historical iterations, and subsequent gradient calculations and direction updates are performed at the prediction point, enabling a Nesterov-accelerated prediction-correction structure. The formula for calculating the prediction point is: ; in, The prediction point corresponding to the nth iteration represents a point on the manifold after incorporating momentum information; For iteration points The manifold retraction operator represents the operational rule for mapping tangent space vectors back to the manifold. The prediction point construction adapts the Nesterov momentum in Euclidean space to the manifold structure, avoiding the problem of iteration points deviating from the feasible region due to direct extrapolation.

[0039] The manifold reduction operator is a smooth mapping that maps tangent space vectors back to the manifold. In this embodiment, it is used to perform manifold mapping operations for momentum prediction and iterative updates, ensuring that the iteration points always satisfy the constraints. In some specific implementations, the manifold reduction operator adopts a normalized reduction form, which can be implemented through simple algebraic operations, has low computational cost, and satisfies the first-order approximation property of manifold mapping. The formula for the normalized reduction operator is: ; in, For point Normalization shrinkage operator at the location; This is the reference point on the manifold, representing the starting point of the retraction operation; Let be a vector in the tangent space, representing the tangent vector that needs to be mapped back to the manifold. The normalization reduction operator's computation process only involves vector addition and norm normalization, which is easy to implement in programming and can guarantee that the output point lies strictly on the unit spherical manifold.

[0040] During the iteration process, a momentum restart mechanism can be introduced. When the inner product of the momentum direction and the current gradient direction is greater than a set threshold, the value of the momentum term is reset to zero, and the value of the auxiliary scalar sequence is reset to its initial value. The momentum direction is the direction after projection of the search direction of the previous iteration into the tangent space, and the current gradient direction is the Riemann gradient direction at the prediction point. The inner product of the momentum direction and the current gradient direction is calculated and compared with a set threshold. When the inner product value is greater than the set threshold, the current momentum accumulation process is terminated, and the momentum update process is restarted; when the inner product value is less than or equal to the set threshold, the current values ​​of the momentum term and the auxiliary scalar sequence are maintained, and the momentum accumulation operation continues.

[0041] Step 3.2. Riemann gradient calculation and processing: First, the Euclidean gradient of the discrete energy function at the prediction point is calculated. Then, the Euclidean gradient is projected onto the tangent space using an orthogonal projection operator to obtain the Riemann gradient. The Riemann gradient represents the gradient information of the energy functional on the manifold and is the basis for constructing the descent direction. In some specific implementations, the Euclidean gradient is the ordinary derivative of the discrete energy function in Euclidean space, reflecting the rate of change of the energy function in each dimension. The Euclidean gradient can be calculated analytically or through numerical difference and is a prerequisite input for generating the Riemann gradient. The formula for calculating the Euclidean gradient is: ; in, Let be the Euclidean gradient at the prediction point of the nth iteration, representing the Euclidean derivative of the discrete energy function at the prediction point; This is the gradient operator, used to compute the gradient vector of a scalar function; It is a discrete energy function; Let be the predicted point for the nth iteration. The Euclidean gradient contains all the information about the changes in the energy function, but it includes the normal component perpendicular to the manifold and cannot be directly used as the descent direction on the manifold.

[0042] In some specific implementations, the Riemann gradient is the result of projecting the Euclidean gradient onto the tangent space, removing the normal component and retaining only the tangential change information. It can be directly used as the steepest descent direction on the manifold. The initial search direction is taken as the opposite direction of the Riemann gradient. The formulas for calculating the Riemann gradient and the initial search direction are: ; ; in, Let be the Riemann gradient at the prediction point of the nth iteration, representing the gradient vector of the discrete energy function on the manifold; This is the Riemann gradient operator, used to compute the gradient of a scalar function on a manifold; For prediction points Orthogonal projection operator at the location; This is the initial search direction, corresponding to the negative Riemann gradient direction. The Riemann gradient naturally lies in the tangent space; moving in the opposite direction of the Riemann gradient can decrease the value of the energy functional, which is the core directional basis for gradient-based optimization methods.

[0043] Step 3.3. Precondition search direction generation process: A positive definite preconditioning operator is constructed, which is an approximation operator of the Riemann-Hessian matrix. This preconditioning operator is used to correct the Riemann gradient, generating a preconditioning search direction. The preconditioning search direction improves the condition number of the optimization problem and accelerates the convergence speed of ill-conditioned problems. In some implementations, the preconditioning Riemann gradient is obtained by combining two tangent space projections with the inverse operation of the preconditioning operator. The preconditioning operator acts on the gradient vector, adjusting the descent step size ratio in each dimension and compressing the spectral distribution range of the objective function. The formula for calculating the preconditioning search direction is: ; in, The precondition search direction for the nth iteration represents the descent direction after precondition correction. For prediction points The preconditioning operator at the given location is used to preprocess the gradient; This is the inverse operator of the preconditioner. The negative sign indicates moving in the opposite direction of the gradient, corresponding to the descent direction of minimizing the energy functional. The two projection operations ensure that the search direction is always within the tangent space.

[0044] The preconditioning operators include the diagonal-dominant operator, the incomplete Cholesky decomposition operator, and the differential principal part approximation operator. The diagonal-dominant operator extracts all elements from the diagonal of the Riemann-Hessian matrix to form a diagonal matrix, using the inverse of this diagonal matrix as the preconditioning operator. The incomplete Cholesky decomposition operator performs a sparse decomposition operation on the Riemann-Hessian matrix, retaining non-zero elements with a set sparsity, generating a lower triangular decomposition matrix, and constructing the preconditioning path based on this lower triangular decomposition matrix. The differential principal part approximation operator extracts the second-order differential principal part operator corresponding to the energy functional and uses the Fourier transform method to achieve fast inversion of the differential principal part operator.

[0045] In some specific implementations, the preconditioner can be constructed in decomposition form, consisting of a combination of Laplace and potential energy terms, approximating the main part of the energy functional Riemann-Hessian matrix. This is suitable for solving energy functionals in quantum physics. The calculation formula for the preconditioner is: ; in, For the Laplace component operator, it corresponds to the preconditioning effect of the kinetic energy part; For potential energy component operators, corresponding to the preconditioning effect of external potential and nonlinear interaction part; is the scaling factor for the Laplace component; is the scaling factor for the potential energy component; is the identity matrix, representing the identity transformation operator; For the Laplace operator, corresponding to the second-order differential operation of the kinetic energy term; Let be the external potential field, representing the spatial distribution function of the external potential energy; The nonlinear interaction coefficient represents the strength of the interaction between particles; The modulus at the prediction point represents the amplitude of the wave function. This decomposition form of the preconditioner can separate the effects of kinetic and potential energy, facilitating the selection of corresponding fast solution methods for different components.

[0046] In some specific implementations, the scaling factor in the preconditioning operator is determined through inner product operations to ensure the adaptability of the operator to the current iteration point, enabling the preconditioning operator to conform to the Hessian properties of the current point. The formula for calculating the scaling factor is: ; The inner product operation is used to calculate the projected component of the operator applied to the current iteration point, ensuring that the scaling factor matches the geometric position of the current manifold point. This form of preconditioning operator can be efficiently inverted using the Fast Fourier Transform without storing the complete matrix, making it suitable for solving large-scale discrete systems.

[0047] Step 4. Line search and iterative update: Step 4.1. Iteration point update processing: A line search is performed along the pre-conditional search direction to determine the iteration step size that satisfies the set conditions. The iteration step size is combined with the pre-conditional search direction to generate an update vector. The update vector is then mapped to a unit spherical manifold using a manifold shrinking operator to generate the next iteration point. Line searches include Armijo line search and Wolfe line search, and the appropriate type can be selected according to the solution requirements. The line search process ensures that the energy functional value decreases sufficiently along the search direction.

[0048] In some specific implementations, the next iteration point is obtained by moving the predicted point along the precondition search direction by a set step size and then back-mapping. Throughout the update process, the iteration point is guaranteed to be located on the manifold surface, requiring no additional normalization correction. The formula for calculating the next iteration point is: ; in, Let be the iteration point of the (n+1)th iteration, representing a point on the updated manifold; For prediction points Manifold shrinkage operator at the location; The iteration step size for the nth iteration is determined by the line search process. After the iterative update operation is completed, the next iteration loop begins until the convergence condition is met, at which point the final solution is output.

[0049] In some specific implementations, line search can use the Barzilai-Borwein step size rule to determine the iteration step size. This rule eliminates the need for multiple function value evaluations, reducing the computational overhead of a single iteration. It includes two different calculation forms, which can be selected according to the actual scenario. The formula for calculating the Barzilai-Borwein step size is as follows: ; in, For the first form of Barzilai-Borwein step size; This is the second form of the Barzilai-Borwein step size; This is the displacement vector from the previous iteration, representing the positional difference between two adjacent iteration points; This is the gradient difference vector from the previous iteration, representing the gradient difference between two adjacent iterations; superscript This represents the vector transpose operation. The two step size forms approximate second-order information from different perspectives, and both can provide reasonable step size values ​​without performing a line search, thus improving the computational efficiency of the iterative process.

[0050] Manifold shrinkage operators include the normalization shrinkage operator and the Riemann exponential mapping operator. The normalization shrinkage operator executes by adding the current iteration point to the tangent space update vector to obtain an intermediate vector, performing Euclidean norm normalization on the intermediate vector, and outputting the normalized vector as the mapping result. The Riemann exponential mapping operator executes by moving a predetermined length along the geodesic direction passing through the current iteration point on the unit spherical manifold, and outputting the final position of the movement as the mapping result. Both the normalization shrinkage operator and the Riemann exponential mapping operator satisfy the first-order approximation property and can map the tangent space vector to the surface of the unit spherical manifold.

[0051] In some specific implementations, the complete solution process is executed using a loop iterative architecture, with each operation module executed sequentially in a fixed order, such as... Figure 2As shown, the solution process begins in the initialization phase, where initial and preceding iteration variables satisfying the unit modulus constraint are given. The relevant parameters of the precondition operator, convergence tolerance, and momentum acceleration parameters are set. After initialization, the iterative loop begins. First, the Nesterov accelerated prediction operation is executed to calculate the momentum direction and generate prediction points. Then, the Euclidean gradient of the energy functional is calculated at the prediction points, and the Riemann gradient is obtained through tangent space projection. Next, the precondition direction is obtained by solving the precondition equation, and the precondition search direction is generated through projection processing. Then, a line search is executed to determine the iteration step size, and a shrinking mapping is used to update the iteration points. The updated iteration points automatically satisfy the unit modulus constraint. After a single iteration, a convergence check is performed, comparing the norm of the Riemann gradient with the set tolerance. If the norm is less than the tolerance, the iteration terminates, and the final solution result is output. If the convergence condition is not met, the iteration number and iteration variables are updated, and the accelerated prediction step is returned to continue the next iteration. This process integrates discrete problem construction, manifold geometry operations, momentum acceleration, preconditioning, and convergence control into a unified closed loop. Each step is connected sequentially, with the output of the previous module serving as the input of the next module, ensuring the orderly execution of the iterative process.

[0052] The preconditioning manifold-accelerated optimization method for spherically constrained energy functionals provided in this embodiment handles spherical constraints through a Riemannian manifold framework, integrating momentum acceleration mechanisms and preconditioning techniques to achieve efficient iterative solutions to constrained optimization problems. This method ensures that the unity modulus constraint is satisfied throughout the iteration process through manifold shrinking operations, eliminating the need for additional projection correction steps, reducing numerical error accumulation, and improving the numerical stability of the iteration process. The introduced Nesterov momentum mechanism utilizes historical iteration information to accelerate convergence, reducing the number of iterations required to achieve convergence accuracy and decreasing the number of evaluations of the energy function and gradient. The introduction of preconditioning operators improves the spectral distribution of the objective functional, reduces the condition number of the optimization problem, and maintains a relatively fast convergence speed even in large-scale discrete systems, mitigating convergence stagnation caused by finer discrete meshes. This method has low single-step iterative computational overhead, does not require storing the complete Hessian matrix, has a small memory footprint, is easily implemented on parallel computing architectures, and is adaptable to high-degree-of-freedom numerical solution scenarios. Meanwhile, the method framework is decoupled from specific physical models; by simply replacing the expressions for the energy functional and its gradient, it can be adapted to different types of spherical constraint optimization problems, exhibiting strong versatility and scalability. Compared with conventional first-order manifold optimization methods, this method has significant advantages in convergence speed and adaptability to ill-conditioned problems, providing a stable and efficient solution path for various energy functional minimization problems with unit modulus constraints.

[0053] The above description is merely a preferred embodiment of the present invention. It should be understood that the present invention is not limited to the forms disclosed herein and should not be construed as excluding other embodiments. It can be used in various other combinations, modifications, and environments, and can be altered within the scope of the concept described herein through the above teachings or related technologies or knowledge. Modifications and variations made by those skilled in the art that do not depart from the spirit and scope of the present invention should be within the protection scope of the appended claims.

Claims

1. A preconditional manifold acceleration optimization method for spherically constrained energy functionals, characterized in that, Includes the following steps: S1. Perform numerical discretization on the target energy functional. The numerical discretization process includes finite difference processing, spectral method processing and finite element processing to generate a discrete energy function and construct a discrete form of spherical constrained optimization problem. S2. For the spherical constraint optimization problem, define the geometric structure of the unit spherical manifold and the corresponding tangent space of the unit spherical manifold. The tangent space is the set of all vectors orthogonal to the current iteration point. Determine the computation rules of the orthogonal projection operator of the tangent space. S3. Introduce an auxiliary scalar sequence to control the momentum coefficient, construct prediction points through manifold shrinkage operators, calculate the Euclidean gradient at the prediction points based on discrete energy functions, project the Euclidean gradient to the tangent space through the tangent space orthogonal projection operator to obtain the Riemann gradient, and modify the Riemann gradient based on the preconditioning operator to generate the preconditioning search direction. S4. Perform a line search along the precondition search direction to determine the iteration step size. Map the update vector to the unit spherical manifold using the manifold shrinking operator to generate the next iteration point. The manifold shrinking operator satisfies the first-order approximation property.

2. The method according to claim 1, characterized in that, The construction of prediction points in step S3 includes the following sub-steps: S31. Update the value of the auxiliary scalar sequence according to the Nesterov acceleration rule. The auxiliary scalar sequence is used to control the degree of momentum decay and generate the momentum coefficient for the corresponding iteration step. S32. The momentum term is generated by linearly combining the momentum coefficient with the displacement vector generated in the previous iteration. S33. The momentum term is mapped from the tangent space to the unit spherical manifold by the manifold shrinking operator, generating the prediction point for the corresponding iteration step.

3. The method according to claim 1, characterized in that, Step S3, calculating the Riemann gradient, includes the following sub-steps: S301. Read the expression of the discrete energy function and calculate the Euclidean gradient of the discrete energy function at the prediction point. The Euclidean gradient is the derivative of the discrete energy function in Euclidean space. S302. Obtain the tangent space orthogonal projection operator corresponding to the prediction point, and project the Euclidean gradient to the tangent space corresponding to the prediction point through the tangent space orthogonal projection operator; S303. Generate the Riemann gradient at the prediction point. The Riemann gradient is the gradient of the target energy functional on the unit spherical manifold.

4. The method according to claim 1, characterized in that, Step S3, which generates the precondition search direction, includes the following sub-steps: S311. Construct a positive definite preconditioning operator. The preconditioning operator is an approximation operator of the Riemann-Hessian matrix. The preconditioning operator satisfies the positive definiteness requirement and has a structure that supports fast inversion calculation. S312. Input the Riemann gradient into the preconditioning operator and perform a linear transformation on the Riemann gradient through the preconditioning operator; S313. Generate a modified precondition search direction, which is the descent direction in the tangent space.

5. The method according to claim 1, characterized in that, Step S4, which generates the next iteration point, includes the following sub-steps: S41. Perform a line search along the precondition search direction. The line search includes Armijo line search and Wolfe line search to determine the iteration step size that satisfies the set conditions. S42. Perform a scalar multiplication operation between the iteration step size and the precondition search direction to generate the update vector for the corresponding iteration step; S43. The update vector is mapped from the tangent space of the current iteration point to the unit spherical manifold using the manifold shrinking operator to generate the next iteration point.

6. The method according to claim 1, characterized in that, In step S1, the spherical constraint optimization problem includes quantum physics computation problem, electronic structure computation problem, numerical algebraic eigenvalue problem, and signal processing constraint optimization problem; the target energy functional corresponding to the spherical constraint optimization problem all contain unit modulus constraints, which require the Euclidean norm of the optimization variable to be equal to a set value; the target energy functional corresponding to the quantum physics computation problem contains kinetic energy term, external potential term, and nonlinear interaction term; The objective energy functional for the electronic structure calculation problem includes a single-electron kinetic energy term, an external potential term, and an electron-electron interaction term; the objective energy functional for the numerical algebraic eigenvalue problem is in the form of a matrix Rayleigh quotient; and the objective energy functional for the signal processing constrained optimization problem is in the form of a signal separation cost function.

7. The method according to claim 2, characterized in that, In step S31, when the inner product of the momentum direction and the current gradient direction is greater than a set threshold, the value of the momentum term is reset to zero, and the value of the auxiliary scalar sequence is reset to the initial value; the momentum direction is the direction after the tangent space projection of the search direction of the previous iteration step, and the current gradient direction is the Riemann gradient direction at the prediction point; the inner product value of the momentum direction and the current gradient direction is calculated, and the inner product value is compared with the set threshold. When the inner product value exceeds the set threshold, the current momentum accumulation process is terminated and the momentum update process is restarted. When the inner product value is less than or equal to the set threshold, maintain the current values ​​of the momentum term and the auxiliary scalar sequence, and continue to perform the momentum accumulation operation.

8. The method according to claim 4, characterized in that, In step S311, the preconditioning operators include a diagonal dominant operator, an incomplete Cholesky decomposition operator, and a differential principal part approximation operator. The diagonal dominant operator extracts all elements on the diagonal of the Riemann-Hessian matrix to form a diagonal matrix, and uses the inverse of the diagonal matrix as the preconditioning operator. The incomplete Cholesky decomposition operator performs a sparse decomposition operation on the Riemann-Hessian matrix, retains non-zero elements with a set sparsity, generates a lower triangular decomposition matrix, and constructs the preconditioning path based on the lower triangular decomposition matrix. The principal part approximation operator extracts the second-order principal part operator corresponding to the energy functional and uses the Fourier transform method to realize the fast inversion calculation of the principal part operator.

9. The method according to claim 5, characterized in that, In step S41, the line search uses the Barzilai-Borwein step size rule to determine the iteration step size, which is calculated using the displacement vector and gradient difference vector from the previous iteration. The Barzilai-Borwein step size rule includes two calculation forms: the first form uses the ratio of the inner product of the displacement vectors to the inner product of the displacement vector and the gradient difference vector as the iteration step size; the second form uses the ratio of the inner product of the displacement vector and the gradient difference vector to the inner product of the gradient difference vector as the iteration step size. The displacement vector is calculated from the position difference between the current iteration point and the previous iteration point, and the gradient difference vector is calculated from the difference between the Riemann gradient of the current iteration point and the Riemann gradient of the previous iteration point. The line search process selects one of the calculation forms to generate the iteration step size.

10. The method according to claim 1, characterized in that, In step S4, the manifold shrinkage operator includes a normalization shrinkage operator and a Riemann exponent mapping operator. The execution flow of the normalization shrinkage operator is as follows: add the current iteration point and the tangent space update vector to obtain an intermediate vector, perform Euclidean norm normalization on the intermediate vector, and output the normalized vector as the mapping result. The execution flow of the Riemann exponential mapping operator is as follows: move a set length along the geodesic direction passing through the current iteration point on the unit spherical manifold, and output the mapping result at the endpoint of the movement. Both the normalization shrinkage operator and the Riemann exponential mapping operator satisfy the first-order approximation property and can map the tangent space vector to the surface of the unit spherical manifold.

Images