Variational formulation for Wannier functions with entangled band structure
Anil Damle Antoine Levitt Lin Lin
physics.comp-ph math.NA physics.chem-ph 65Z05, 82D25, 65F30, 65K10


Effective Building Block Design for 
Deep Convolutional Neural Networks using Search
    Jayanta K Dutta1, Jiayi Liu1,
	Unmesh Kurup1,
	Mohak Shah2

	1Bosch Center for Artificial Intelligence, Sunnyvale, CA, USA

	2University of Illinois at Chicago, Chicago, Illinois, USA 

	

    December 30, 2023
================================================================================================================================================================================================




  Wannier functions provide a localized representation of spectral
  subspaces of periodic Hamiltonians, and play an important role for
  interpreting and accelerating Hartree-Fock and Kohn-Sham density
  functional theory calculations in quantum physics and chemistry. For
  systems with isolated band structure, the existence of exponentially
  localized Wannier functions and numerical algorithms for finding
  them are well studied. In contrast, for systems with entangled band
  structure, Wannier functions must be generalized to span a subspace
  larger than the spectral subspace of interest to achieve favorable spatial
  locality. In this setting, little is known about the theoretical
  properties of these Wannier functions, and few algorithms can find
  them robustly. We develop a variational formulation to compute these
  generalized maximally localized Wannier functions. When
  paired with an initial guess based on the selected columns of the
  density matrix (SCDM) method, our method can robustly find Wannier
  functions for systems with entangled band structure. We formulate
  the problem as a constrained nonlinear optimization problem, and
  show how the widely used disentanglement procedure can be
  interpreted as a splitting method to approximately solve this
  problem. We demonstrate the performance of our method using real materials including
  silicon, copper, and aluminum. To examine more precisely the
  localization properties of Wannier functions, we study the free
  electron gas in one and two dimensions, where we show that the
  maximally-localized Wannier functions only decay algebraically. We
  also explain using a one dimensional example how to modify them to obtain
  super-algebraic decay.





Wannier function, Localization, Entangled band, Metallic system,
Variational method, Optimization, Free electron gas


myheadings
plain



§ INTRODUCTION


Localized representations of electronic wavefunctions have a wide range
of applications in quantum physics, chemistry, and materials science.
In an effective single particle theory such as Hartree-Fock theory
and Kohn-Sham density functional theory
(KSDFT) <cit.>, the electronic
wavefunctions are characterized by eigenfunctions of single particle
Hamiltonian operators. These eigenfunctions generally have significant
magnitude in large portions of the computational domain. However, the physically meaningful quantity is
not each individual eigenfunction, but the subspace spanned by the
collection of a set of eigenfunctions. This is often referred to as the Kohn-Sham subspace, and it is often possible to reduce the computational complexity of various methods by using an alternative, localized representation of the subspace.


Wannier functions provide one such localized representation of the Kohn-Sham
subspace. They require
significantly less memory to store, and are the foundation of
so-called “linear scaling
methods” <cit.>
for solving quantum problems.  They can also be used to analyze chemical
bonding in complex materials, interpolate the band structure of
crystals, accelerate ground and excited state electronic structure
calculations, and form reduced order models for strongly correlated
many body systems <cit.>.


Wannier functions are not uniquely determined, and depend on a
choice of gauge (a rotation among the occupied states), which strongly
influences their localization. For periodic systems with an isolated
band structure, the localization properties of Wannier functions have
been studied extensively
<cit.>.
Interestingly, the existence of localized Wannier functions in this
case is characterized by a topological invariant. For physical systems
without magnetic field (invariant under time-reversal symmetry), this
topological invariant is trivial, and it is known that there exists a
gauge leading to Wannier functions with exponential decay
<cit.>. In this
setting, efficient numerical algorithms have been developed to compute
these exponentially localized functions
<cit.>.
In particular, the widely used maximally-localized Wannier function
(MLWF) procedure minimizes the variance (or
“spread”) <cit.> over all
possible choices of gauge to obtain localized Wannier functions. In
the insulating case, it is known that minimizers of this spread are
exponentially localized <cit.>.

The situation becomes significantly more challenging for systems with
entangled band structure. Entangled band structure arises in metallic
systems, but also in insulating systems when conduction bands or a
selected range of bands are to be localized. A straightforward
definition of Wannier functions requires the set of all Wannier
functions to exactly span the selected spectral subspace. However,
such Wannier functions are known to decay slowly in real space.
Therefore, the definition of Wannier functions has been generalized to
refer to functions spanning a subspace larger than, but containing the given
entangled spectral subspace, referred to as a “frozen window”
<cit.>. This is useful for instance in band
interpolation, where the additional Wannier functions give rise to
extra bands that can simply be ignored. Finding such generalized Wannier
functions numerically is considerably more complex, and few algorithms
in the literature accomplish this task in a robust fashion.
Furthermore, little is known theoretically about the localization
properties of the constructed generalized Wannier functions. In order
to be consistent with the terminology in the physics literature, we
will refer to these generalized Wannier functions simply as Wannier
functions, unless otherwise noted.



In this paper, we develop a variational formulation for finding Wannier
functions in the entangled setting. We formulate the problem as a nonlinear constrained
optimization problem. Practical Wannier function calculations indicate
that such nonlinear optimization problems can have many local minima.
Hence the solution can strongly depend on the initial guess, and the difficulty of constructing a good initial guess is often a significant impediment to
finding Wannier functions in a robust fashion. In order to avoid
being trapped at undesirable local minima, we use the recently
developed selected columns of the density matrix (SCDM) methodology to
construct the initial guess for our variational formulation. This
strategy is applicable to both the isolated
case <cit.> and the entangled
case <cit.>.


Our variational formulation can be obtained in several theoretically
equivalent constructions.  We find that one of these formulations yields
the so-called partly occupied Wannier
functions <cit.>, and can be solved efficiently using
standard numerical algorithms for minimization under orthogonality constraints.
Our formulation also reveals that the widely used “disentanglement”
procedure <cit.> can be viewed as a splitting
method for solving the constrained optimization problem, which only
performs a single alternation step between the two pieces of the
objective function, and therefore does not achieve a global minimum of
the spread. We verify the performance of the variational
formulation with real materials such as silicon, copper, and aluminum.
In these examples, we find that the fully converged variational
formulation consistently provides orbitals with a smaller spread than
that from the disentanglement procedure, and is more robust to the
choice of initial guess. We also find that the difference between the
variational formulation and the disentanglement procedure is often
small when used for band structure interpolation.





The variational formulation allows us to study the decay properties of
Wannier functions for metallic systems. We present the localization
properties of generalized Wannier functions for the free electron gas
in one and two dimensions. We find that minimizers of the spread
exhibit a weak algebraic decay, related to singularities that we
identify in  space. This slow decay is shown to be not a
fundamental property of disentangled Wannier functions, but rather a
consequence of the fact that minimizing the spread only imposes finite
second moments (or square-integrable first derivatives in 
space). In particular we show in one dimension how to modify the
maximally-localized Wannier functions to obtain super-algebraic decay.

The rest of the paper is organized as follows.  We first introduce several background topics such as Bloch-Floquet theory, Wannier functions,
and the SCDM methodology in
section <ref>. We then present a variational formulation for
Wannier functions in section <ref>, and discuss the
relation between our variational formulation and existing methods.
Numerical results for real materials and for the free electron gas are
given in sections <ref> and  <ref>, followed by conclusion and discussion
in section <ref>. Some of the technical details related
to the implementation of the variational formulation are given in the Appendix <ref>.



§ PRELIMINARIES




 §.§ Bloch-Floquet theory


We first briefly review
Bloch-Floquet theory for crystal structures. 
The Bravais lattice with lattice vectors _1,_2,_3∈^3 is defined as

    𝕃 = { = ∑_i=1^3 n_ia_i,    n_1,n_2,n_3∈ℤ},

and the lattice vectors define a unit cell in the Bravais lattice

    Γ ={ = ∑_i=1^3 c_ia_i| -1/2 ≤ c_1,c_2,c_3< 1/2}.

The Bravais lattice induces a reciprocal lattice denoted
𝕃^*, which is the support of the Fourier transform of
𝕃-periodic functions. The lattice vectors of 𝕃^*
are denoted by
_1,_2,_3, with _i·_j = 2πδ_ij. A unit cell of the
reciprocal lattice is selected and called (with some abuse of
language) the Brillouin zone, and is defined as

    Γ^* ={ = ∑_i=1^3 c_ib_i| -1/2 ≤ c_1,c_2,c_3< 1/2}.


For a potential V that is real-valued and 𝕃-periodic,   

    V( +) = V(),   ∀∈^3, ∈𝕃,
 
we consider the Schrödinger operator in ℝ^3

    H = -1/2Δ + V.



The Bloch-Floquet theory allows us to relabel the spectrum of H using two indices (n,), where n∈ℕ is the band
index, and ∈Γ^* is the Brillouin zone index. The
generalized (not square-integrable) eigenfunction ψ_n,() is known as a Bloch orbital and 
satisfies

    Hψ_n,()=ε_n,ψ_n,().

Importantly, ψ_n, can be decomposed as 

    ψ_n,() = e^· u_n,(),

where u_n,() is a periodic function with respect to
𝕃. Eigenpairs (ε_n,,u_n,) can therefore be 
obtained by solving the eigenvalue problem

    H() u_n, = ε_n, u_n,,    n ∈ℕ,   ∈Γ^*,

where H() = 1/2 (-∇+ )^2 + V().
For each , the eigenvalues ε_n, are ordered
non-decreasingly, and {ε_n,} as a
function of  for a fixed n is called a
band. The set of all eigenvalues is called the
band structure of the crystal and characterizes the spectrum
of the operator H. 




If ε_N+1, > ε_N, for all ∈Γ^*,
then the first N bands are isolated. This is for instance the case
in the occupied bands of an insulator. 
When the gap condition ε_N+1, > ε_N,
does not hold, the band
structure becomes entangled. Entangled band structure appears not only
in metallic systems, but also insulating systems when a Wannier
representation of part of the conduction bands is required.















 §.§ Wannier functions



For simplicity, we first consider systems with isolated first N bands—an assumption we will drop towards the end of this section. 
Rotating the set of functions {ψ_n,} by an arbitrary
unitary matrix U()∈^N× N, we can
define a new set of functions

    ψ_n,() = ∑_m=1^Nψ_m,() U_m,n(),   ∈Γ^*.

A given set of of such matrices {U()}_∈Γ^* is called a gauge.

For each ,
we consider the density matrix P(), which is the projector on the
the eigenspace corresponding to the first N eigenvalues of H()

    P() = ∑_n=1^N|ψ_n,⟩⟨ψ_n,| =∑_n=1^N|ψ_n,⟩⟨ψ_n,|.

Importantly, for each , the density matrix P() is
gauge-invariant. If 𝒞 is a contour in the complex plane enclosing the eigenvalues
ε_1,,…, ε_N, (and only those), then the Cauchy
integral formula yields an alternative representation of P()

    P() = 1/2π∫_𝒞1/λ - H()λ.

Since H() is analytic, it
follows that so is P().

Given a choice of gauge, the Wannier functions are defined as <cit.>

    w_n,() = 1/|Γ^*|∫_Γ^*ψ_n,()
      e^-·,   ∈^3, ∈𝕃,

where |Γ^*| is the volume of the first Brillouin zone. This represents a unitary transformation from the family
(ψ_n,)_n=1,…,N,∈Γ^* to
(w_n,)_n=1,…,N,∈𝕃. In particular, the
Wannier functions w_n, are orthogonal to each other and span
the same space as the range of the total density matrix
1/|Γ^*|∫_Γ^* P(). They are also
translation invariant: w_n,() = w_n(-).

For insulating systems, in the absence of topological obstructions,
there exists a gauge such that ψ_n, is analytic and
𝕃^*-periodic in , implying that the Fourier
transform of w_n, is analytic, and therefore that each
Wannier function decays exponentially as
→∞ <cit.>.
The Wannier localization problem is reduced to the problem
of finding a gauge {U()} such that w_n,0 is localized, or, equivalently, that ψ_n, is smooth with respect
to . This can be done by minimizing the “spread
functional” <cit.>

    Ω[{U()}] = ∑_n=1^N∫w_n,0()^2^2  - | 
       ∫w_n,0()^2|^2.

Here w_n,0 depends on U() through ψ_n, as
in Eq. (<ref>).  This problem is usually
solved by a minimization algorithm such as steepest descent or
conjugate gradient with projections at each step to respect the constraints
that U() must be unitary <cit.>.









For systems with entangled band structure, the density matrix P()
as defined in Eq. (<ref>) is no longer smooth with
respect to . As a result, there is no choice of gauge U()
that leads to a set of rotated Bloch orbitals that is smooth with
respect to , and Wannier functions defined strictly according to
Eq. (<ref>) will then decay very slowly in real
space <cit.>. In order to enhance the localization
properties of Wannier functions, the definition of Wannier functions
has been generalized so that the spectral subspace interest is only a
proper subspace of that spanned by Wannier
functions <cit.>. In the physics literature,
the spectral subspace is described by a “frozen window” along the
energy spectrum, and the Wannier functions are linear combinations of
orbitals from a larger set described by an “outer window”.








More specifically, we first fix a number of bands N_o that
determines the outer
window[To simplify the exposition we assume a constant number of bands in the outer window, but this can be relaxed to a variable number of bands N_o().] and then proceed to look for Wannier functions
built out of ψ_1,,…,ψ_N_o,. Next, for each  point we fix a
set of frozen bands 𝒩_f() ⊂[N_o], and let N_f() = |𝒩_f()|. 𝒩_f() are often defined as the bands within a fixed energy window that we
will try to reproduce. Correspondingly, we
define the frozen density matrix as 

    P_f() = ∑_n∈𝒩_f()|ψ_n,⟩⟨ψ_n,|,

which is the projection onto the states within the frozen energy
window. Again, the frozen density
matrix as defined in Eq. (<ref>) is not smooth with
respect to .








We now seek to construct a set of N_w Wannier functions
that span the subspace defined by the range of
1/|Γ^*|∫_Γ^* P_f(). We introduce
the gauge matrices U()∈^N_o× N_w with orthonormal
columns, with | N_f()|≤ N_w≤ N_o, such that













    ψ_n,() = ∑_m=1^N_oψ_m,() U_m,n(),   ∈Γ^*,  n=1,…,N_w.

This may equivalently expressed in matrix form as Ψ() = Ψ() U().
This choice of gauge also induces a density matrix of rank N_w for each  defined as

    P_w() = ∑_n=1^N_w|ψ_n,⟩⟨ψ_n,| = Ψ()
      Ψ^*() = Ψ() U() U^*() Ψ^*().

Note that unlike the case with isolated band structure where N_o=N_w, here N_w≠ N_o implies that U()
U^*() ≠ I_N_o.  Furthermore, since the set of orbitals in the frozen window is only a subset of all
possible orbitals, in general the projectors P_w and P_f do
not span the same space. In
order to ensure that our Wannier functions span the same subspace as the
subspace associated with the frozen window, we require that

    P_w() P_f() = P_f(),   ∀∈Γ^*.












§ VARIATIONAL FORMULATION FOR WANNIER FUNCTIONS WITH ENTANGLED
BAND STRUCTURE


We now proceed to develop a variational formulation for Wannier
functions. First, we illustrate how to encode the desired constraints
when paired with the aforementioned spread functional. Subsequently, to facilitate numerical solution of the optimization problem, we refine how the constraints are expressed. Lastly, we discuss the relation the existing disentanglement procedure to our formulation and discuss how we construct an initial guess using the SCDM methodology.



 §.§ Formulating the optimization problem

Without loss of generality, for the following discussion we assume that the
frozen orbitals (Ψ_f) are always ordered before the rest of the
orbitals (Ψ_r). In terms of the notation from the previous
section, this means that for each  the frozen orbitals are
represented by the set 𝒩_f() = {1,2,…,N_f()} with
N_f() simply representing the number of frozen orbitals per
-point. Now, we may partition the orbitals and the gauge using
the following block form

    Ψ()  = [ Ψ_f() Ψ_r() ],
          U() =
      [ U_f(); U_r() ].



The matrices U_f() and U_r() are of size
N_f() × N_w and (N_o - N_f()) × N_w, 
encoding the weight assigned to the frozen
subspace and its complement in the Wannier functions, respectively. The condition that the Wannier functions represent the
frozen bands as in Eq. (<ref>) can conveniently
be expressed in terms of these matrices as follows.


  The following statements are equivalent:
  
    
  * P_w() P_f() = P_f().
    
  * U_f() U^*_f() = I_N_f().
    
  * U_f() U^*_r() = 0 and U_f() has full row
      rank.
    
  * U()=[ I_N_f()       0;       0     Y() ] X(),  where X() is a unitary matrix of size N_w× N_w,
  and  Y() is a matrix with orthogonal columns of size
  (N_o-N_f()) × (N_w-N_f()).
  


  Since each  point is treated independently, for simplicity we
  drop the  dependence in the proof below.

   1 ⇔ 2: From the definition of
   P_w,P_f we have
   
    P_wP_f   = Ψ[ U_f U_f^* U_f U_r^*;  U_rU_f^*  U_rU_r^* ][ I 0; 0 0 ]Ψ^* = Ψ[ U_f U_f^*         0;  U_rU_f^*         0 ]Ψ^*
    
         P_f   = Ψ[ I_N_f     0;     0     0 ]Ψ^*

   The result follows since Ψ has orthogonal columns.


   2 ⇔ 3:
     From the partition of U we have
     
    U U^* =
           [ U_f U_f^* U_f U_r^*;  U_rU_f^*  U_rU_r^* ]

   If U_f U_f^* = I_N_f, since UU^* is a projector, it
   follows that U_rU_f^* = 0. On the other hand, if
   U_rU_f^* = U_fU_r^* = 0, then the fact that UU^* is a
   projector implies that U_fU_f^* is a projector as well. Since
   it has full row rank, it must therefore be the identity matrix.

   3 ⇔ 4: 
   If 3 is true, then
   
    P_U = U U^* = 
         [    I_N_f        0;        0 U_rU_r^* ].

   Then U_rU_r^* is a projector with rank N_w-N_f, or equivalently
   
    U_rU_r^* = Y Y^*

   for some (N_o-N_f) × (N_w-N_f) matrix Y with
   orthogonal columns. Then
   
    U = P_U U = 
         [ I_N_f     0;     0 Y Y^* ] U = 
         [ I_N_f     0;     0     Y ] X,

   where 
   
    X=[ I_N_f     0;     0   Y^* ] U,

   and it can be readily verified that X is unitary. The reverse
   direction is obvious.



Proposition <ref> gives us
various concise ways to impose the desired condition on the span of
Wannier functions, and we may for instance consider using condition 2. Since the smoothness requirement for ψ_n, with respect to the
Brillouin zone index  can be realized by minimizing the spread
functional (<ref>), finding the desired smooth gauge U()
can be recast as the following constrained optimization problem:


    inf_{U()}     Ω[{U()}]
    s.t.      U^*() U() = I_N_w,   
      U_f() U_f^*() = I_N_f().


The difficulty at this stage is that numerical optimization of (<ref>) with
respect to these constraints may not be easy. In particular, the set of
matrices U satisfying U^*U = I_N_w and U_fU_f^* = I_N_f does not 
necessarily possess a smooth manifold structure. This complicates the application
of standard methods for the minimization of functions over
orthogonality constraints. 

On the other hand, the condition 4 in
Proposition <ref> represents U in a factorized form,
hereinafter referred to as the (X,Y) representation.
This representation of the matrix U() gives rise to Wannier
functions composed of the N_f() functions in the
frozen window, and another set of N_w - N_f()
functions, encoded by the matrix Y(). This Y encapsulates all
the necessary information about the projector P_w = ΨΨ^*. The unitary
X() matrix mixes these N_w Wannier functions amongst themselves to
produce a smooth gauge. In the (X,Y) representation,
the variational formulation (<ref>) can be written
as

    inf_{X(),Y()}     Ω[{U()}],
    s.t.      U() = [ I_N_f()       0;       0     Y() ] X(), 
          X^*() X() = I_N_w,
          Y^*() Y() =
      I_N_w-N_f().

This optimization problem is equivalent to the “partly occupied Wannier
functions” <cit.>. This also directly
generalizes the maximally localized Wannier functions
procedure <cit.> by Marzari and Vanderbilt for the
isolated case.



The (X,Y) representation is a redundant representation, and a given U
can be reproduced by many pairs (X,Y). However,
in contrast to the formulation (<ref>), the constraint in
Eq. (<ref>) defines a Riemannian manifold where
X() and Y() are independent matrices with orthogonality
constraints. This allows us to use standard algorithms for the
minimization of differentiable functions on Riemannian manifolds to
solve the problem. We refer to Appendix <ref> for the details
of the computation of the gradient of the objective function Ω.



 §.§ Implementation


For our implementation, we modified the Julia <cit.>
library  for unconstrained optimization to
accommodate constraints represented by Riemannian manifolds
<cit.>. Our modifications
have been integrated into that library and are available online
[<https://github.com/JuliaNLSolvers/Optim.jl>]. For the
numerical tests that follow we used the limited-memory BFGS algorithm
<cit.> with Hager-Zhang line search
<cit.>, which gave the best performance compared to other
readily-available algorithms (steepest descent, conjugate gradient,
BFGS) and line searches (fixed step, backtracking).


In order to generate the initial guess for numerical optimization, we
need to convert a given matrix U to a pair (X,Y) that
parametrizes it. It will also be useful to consider matrices U that
do not satisfy the constraints U^*U=I_N_w and U_fU_f^*=I_N_f exactly but only approximately. 


This will allow us to project U to the admissible set that satisfy
these constraints.

To find a pair (X,Y) that represents a given U, we
first choose Y to minimize the error on the projector UU^*
measured by the Frobenius norm via

    inf_Y^*Y=I‖ UU^*-[ I_N_f     0;     0 Y Y^* ]‖_F^2.

A solution to this problem can be computed using the eigenvalue decomposition

    U_rU_r^* = V S V^*.

When U satisfies the constraints, U_rU_r^* is a
projector of rank N_w-N_f, and we can choose Y to be the
columns of V corresponding to the N_w-N_f non-zero eigenvalues
of S. When U does not satisfy the constraints, we pick Y as the
columns of V corresponding to the largest N_w-N_f eigenvalues
of S.

Once Y is computed, we can find X that minimizes the error on U:

    inf_X^*X=I‖ U-[ I_N_f     0;     0     Y ] X‖_F^2.

When U satisfies the constraints, the solution of this problem
is simply

    X
    =  [ I_N_f     0;     0   Y^* ] U.
 
Otherwise, the solution can again be obtained via the singular value decomposition

    [ I_N_f     0;     0   Y^* ] U = V_ lSV_ r^*,

and setting X = V_ lV_ r^*. This step is also
called the Löwdin orthogonalization procedure <cit.>.



 §.§ Relation to disentanglement

Our variational formulation also gives us a concise way to understand the
“disentanglement” procedure of Souza, Marzari and
Vanderbilt <cit.>, in which the spread
functional is split into two parts

    Ω[{U()}] = Ω_I[{U()}] +
      Ω[{U()}].

Here Ω_I is called the gauge-invariant part (depending on P_w,
and hence only on Y Y^*), and Ω is called the gauge
dependent part (depending on X).
Instead of optimizing (<ref>)
directly, <cit.> proposes to use a two step
procedure. First one optimizes the gauge-invariant part only:

    inf_{Y()}  Ω_I[{U()}].

This is numerically expedient as Ω_I only depends on
Y Y^*. In fact, it is analogous to minimization problems in electronic
structure (for instance the Hartree-Fock model), where one minimizes
the energy, which only depends on the spectral subspace, over all
possible orthogonal orbitals. The authors in
<cit.> accordingly obtain a nonlinear eigenvalue
problem as the first-order optimality conditions, which they solve
using a damped self-consistent field (SCF) iteration.

Once {Y()} is obtained, it is fixed and so is the projector
P_w. A second minimization problem

    inf_{X()}  Ω[{U()}]

is then solved with respect to the gauge matrix X(). This
optimization problem is of
the same nature as the one for an isolated set of bands.

The total spread from the above two-step
procedure is necessarily larger or equal to the global minimum
of (<ref>). Interestingly, although the optimal spread can be
substantially lower than that found by the two-step disentanglement
procedure, numerical experiments show that the quality of Wannier
interpolation, measured for instance by the qualify of band structure
interpolation, is often similar in both cases.



 §.§ Selected column of the density matrix


While the primary purpose of this work is to introduce and analyze a
variational formulation of Wannier functions, both the objective
function and the constraints are nonlinear, and hence there may exist
multiple local minima.  It is practically important to seed
such methods with a good initial guess. Here, we summarize the
recently developed unified methodology for Wannier localization of
entangled band structure <cit.> based on the selected
columns of the density matrix (SCDM) methodology <cit.>. Importantly, this method is direct and robust—no initial guess is required and it will generate valid output—and thus may be reliably used to generate an initial guess.

The SCDM method for entangled band structure first constructs a
quasi-density matrix

    f(H()) = ∑_n=1^N_of(ε_n,)
      |ψ_n,⟩⟨ψ_n,|,
  
For insulating systems f would be 1 on the occupied bands and
0
otherwise, yielding the projector P as before. For entangled band
structure however, the function f(·) is chosen to be large on the bands of
interest and decays rapidly, but smoothly, away from them
<cit.>. The SCDM algorithm constructs a gauge
by selecting a common set of columns of the -dependent
(quasi-)density matrix f(H()). In practice, it is often sufficient to
select these columns based on an “anchor” point denoted _0—generically chosen to be the so-called Gamma-point (0,0,0)^T.

We now briefly outline the SCDM method and refer the reader to
<cit.> for more details. Let
Ψ_∈ℂ^N_g× N_o be the matrix with
orthogonal columns that represents {ψ_n,()} on
a discrete grid in the unit cell, and E() =
diag[{ε_n,}_n=1^N_o]
be a diagonal matrix encoding the corresponding eigenvalues.
SCDM identifies N_w columns of f(H()) based on the leading N_w
columns of the permutation matrix Π, computed via the
QR factorization with column pivoting (QRCP) procedure

    Ψ__0^*Π = QR.

This set of columns is denoted by C =
{_n}_n=1^N_w⊂Γ. Now, for each -point define Ξ()∈ℂ^N_o× N_w as

    Ξ_n,n'() = f(ε_n,)ψ_n,^*(_n').
 
It is expected that 

    ψ_n,() = ∑_mψ_m,()Ξ_m,n()

is smooth with respect to . Therefore, if the singular values of
Ξ() are uniformly bounded away from 0 in the Brillouin zone,
U() constructed via Löwdin orthogonalization <cit.> of Ξ() has
orthogonal columns, and yields {ψ̃_n,}
that are smooth with respect to . 

In this framework the frozen bands are
not represented exactly, so prior to use in our optimization procedure we must convert from U to a pair (X,Y) via our aforementioned scheme. Since the projector on the
frozen set varies discontinuously with , this procedure does not
produce a continuous gauge. However, if the function f is chosen appropriately, it should be close to one. This is further corroborated by the quality of band interpolation, despite the substantially larger spread, we achieve in the numerical results section using the SCDM initial guess without explicitly freezing any bands.







§ REAL MATERIALS

We first consider the performance of our variational formulation for
several real materials. This includes valence and conduction bands of
silicon (semiconductor), conduction bands of copper (metal), and valence
bands of aluminum (metal). We always start with the aforementioned SCDM
based initial guess.  We compare the result obtained from the
variational formulation to that obtained from the disentanglement
procedure in , as well as the result obtained directly from the
SCDM initial guess without further refinement. 

In these experiments, the choice of the
parameters of the SCDM procedure are chosen to yield good baseline band
structure interpolation. However, they are not “optimized” to minimize band structure interpolation error. These experiments often show how the two optimization methods are
comparable, though in some situations we are able to find Wannier functions with
smaller spreads using our variational method even given the same initial guess. One interesting
point that we will see play out throughout our examples is that the
value of the spread and band interpolation quality may not be directly
connected, i.e. 
Wannier functions with considerably different values of spread can yield
qualitatively similar interpolation error.

All of the -point calculations and reference band structure
calculations were performed with
 <cit.>. The SCDM initial guess was
constructed using the code available
online[<https://github.com/asdamle/SCDM>]. Our new variational
formulation was implemented in the Julia language and is available
online[<https://github.com/antoine-levitt/wannier>].




 §.§ Silicon

Here, we compute the lowest 16 bands of Silicon on an 8× 8× 8 -point grid and then proceed to
compute eight Wannier functions. This includes the four valence bands
and four additional low lying conduction bands. For the SCDM procedure
we use μ = 11.0 and σ = 2.0 with f corresponding to
“entangled case 1” in <cit.>—a complementary error function. For and our
method,
we freeze bands below 12 eV and set the outer window
maximum at infinity. For the prescribed convergence criteria
of 1× 10^-10 on the spread was reached after 225
disentanglement iterations and 95 spread reduction iterations, and for
our method after 149 iterations.

Figure <ref> shows the band structure interpolation
using the three methods. We see that for the four
valence bands all three methods perform very well. Furthermore, while
there are differences in the interpolation of the conduction bands, no
one method clearly outperforms the others. As expected, if we consider
the total spread (see Table <ref>) of the final localized
orbitals, our variational formulation yields the most compact orbitals.

In Table <ref> we report the per orbital spread for
each method, and observe that to the number of significant digits
reported all the orbitals found by our method have the same spread (they
do not vary until the fourth decimal place). In contrast, seems
to converge to two distinct sets of orbitals with slightly different
spreads. In Figure <ref>, we illustrate the differences
between the orbitals found with our variational method and those found
via .  We observe that the orbitals obtained from our
variational formulation resemble more closely to sp^3 hybridized
orbitals centering around each Si atom, as indicated from 
chemical intuition.  Interestingly, by using the output of our variational
method as input to we are able to force to converge to
the same point as our method. Unfortunately, it is difficult to pinpoint
a specific cause for the apparent convergence of to a worse
local minimum in this setting.
































 §.§ Aluminum

We now repeat the same experiments as before, but for the valence bands
of the aluminum system with an 8× 8× 8 -point grid. Specifically, we start with six bands and seek
four Wannier functions. For the SCDM procedure we use μ = 8.42 and
σ = 4.0 with f corresponding to “entangled case 1” in
<cit.>—a complementary error function. Once again, was run until convergence at 10^-10 for both the disentanglement and spread minimization, and reached that threshold after 2,437 and 91 iterations. Our method converged to a spread reduction tolerance of 10^-10 after 138 iterations. As we observe in Figure <ref>, all three methods once again perform well—particularly below
the Fermi energy. The final spreads for each of the three methods are
reported in Table <ref> along with the spreads of each orbital. While both and our
method improve upon the spread of the SCDM initial guess, we do find
slightly smaller spread with our optimization procedure. Here, bands
below 11.6 eV were frozen in both and our method, and the outer
window was set to ∞.


















We also consider the convergence of the band interpolation error in this setting, looking at both maximum error and RMSE. As before, for both our method and we freeze bands below 11.6 eV and set the outer window to infinity. For all three methods we then measure the error of band interpolation at or below the Fermi energy (8.42 eV). In all cases, we capped at 5,000 disentanglement and spread reduction iterations and considered it converged at a tolerance of 10^-10. Similarly, we considered our method converged the objective function changed by less than 10^-10 between successive iterates, or we reached 5,000 iterations. Typically our method took 100-250 iterations to converge, whereas the disentanglement would hit the iteration cap and spread reduction took 90-400 iterations to converge. Figure <ref> shows broadly similar behavior for all three methods, though generally the two optimization methods do noticeably improve on the SCDM initial guess as more -points are used. We expect that asymptotically the optimization based methods should perform better. However, given the relatively small number of grid points per direction and the complexity of the band structure, it seems we are still in the pre-asymptotic regime.







 §.§ Copper

Lastly, we illustrate the behavior of our method when used to
interpolate seven conduction bands of copper around the Fermi energy.
Because several bands that we do not wish to track cross through the
energy window, it is not easy to simply look at band interpolation
error. Rather, we place emphasis on the qualitative behavior of the
interpolation and the orbital spread yielded by our variational
formulation. Interestingly, in this case we observed a high
sensitivity of to the SCDM initial guess based on the
parameters (corresponding to “entangled case 2” in
<cit.>—a Gaussian) and a contrasting
robustness of our variational method. In all cases we used a frozen
window of 13.5 to 17 eV and no outer window for both and our variational method.

When fixing the parameter μ = 15.5 and varying σ from 3.0
to 6.0, both SCDM and our variational method robustly generated good
band interpolation. However, for a range of σ with
disentanglement failed to reach convergence. To sweep over several
values of σ we limit to 5,000 iterations each for the
disentanglement and spread minimization, and we limit our variational
method to 1,000 iterations. 


For σ = 5.0 where we
observed particularly bad performance of (see below), we let it run for
100,000 iterations of each the disentanglement and spread. While the
disentanglement procedure converged after roughly 20,000 iterations,
the spread minimization failed to converge even after 100,000
iterations. While this does not guarantee that the local minimum
may eventually find is poor and could simply be the
optimization algorithm behaving poorly, we feel this is a reasonable
comparison to make even without determining that it has
converged. Interestingly, passing the output of in this setting to our variational method we were able to converge to a good gauge.




Figures <ref> and <ref> show
the band interpolation of the three methods in the case where σ =
4.0 and σ = 5.0. We also report the individual spreads of the orbitals in Tables <ref> and <ref>.  We observe that in both cases the SCDM initial
guess and our optimized solution yield good band interpolation within
the frozen window. In
contrast, does not find a good local optima in the latter case
and this results in poor interpolation quality. We further investigate
this behavior in Table <ref>, where we report the total
final spread of the methods as we vary σ. We see that for two of
the parameter values failed to find a local optimum that is close
to what our variational method finds.


















§ FREE ELECTRON GAS

We now investigate the decay property of the generalized Wannier
functions in real space. This is hard to
investigate numerically for real materials because of the very large
number of  points needed to see the asymptotic decay, as
evidenced in the previous section. For
instance, in <cit.> the authors report a
fast (consistent with exponential) decay up to grid sizes of
15 × 15 × 15, although the convergence appears to slow
down after that. Another problem is that it is hard to find realistic
systems in one or two dimensions on which disentangled Wannier systems
make sense. The free electron gas,  V = 0, is explicitly solvable
and poses a very interesting benchmark for disentanglement, even in
one or two dimensions. In this section, we therefore apply the
methodology of Section <ref> to d-dimensional free
electron gas. We simply define the lattice as
𝕃 = 2πℤ^d, so that
𝕃^* = ℤ^d, and we let the Brillouin zone be the
set [0,1)^d.

In this case the eigenfunctions of the operator H() = (-∇
+ )^2 are given by

    v_()    = 1/√(|Γ|) e^·

for  in the reciprocal lattice ℤ^d. The corresponding
eigenvalues are ε_, = | + |^2. Notice that
since ε_, is the squared distance from  to
-, the dispersion relation is the set of squared distances of
 to the points of the reciprocal lattice ℤ^d.

We order the eigenfunctions by ordering the eigenvalues in a
non-decreasing order. We let
u_n, = e__n, where ε__n, is the
n-th eigenvalue of H(), this choice being arbitrary in the
presence of degeneracies. The matrix elements
M_mn^, = ⟨ u_m,, u_n,+⟩ of
overlap between neighboring  points used in
the optimization process <cit.> then assume a
particularly simple expression: M_mn^, = 1 if u_m
and u_n,+ are associated with the same , and M_mn^, = 0
otherwise. In particular, this matrix differs from the identity (it is a permutation matrix) near eigenvalue crossings, where
ε_, = ε_', with ' ≠.

The free electron gas also makes it particularly easy to compute the
Wannier functions through their Fourier transforms. For a Wannier
function given in -space by

    ψ_n() = ∑_n ∈ e^·
     v__m()U_mn(),

it holds that

    w_n,0() = ∫_[0,1]^dψ_n()  = 1/(2π)^d∫_[0,1]^d∑_m ∈ U_mn()e^ (+_m)·,

from which it follows that U_mn() is simply the Fourier transform
of w_n,0 at frequency ξ = _m+.

Note that the free electron gas possesses a large number of symmetries. As a
consequence, it has a number of properties that are not expected from
generic systems. For instance, eigenvalue crossings are numerous (of
co-dimension 1,  points in 1D, lines in 2D and planes in 3D), while
they are expected from the von Neumann-Wigner theorem
<cit.> to be rare. This theorem predicts that, in
the absence of particular symmetries, the crossing of eigenvalues of a
Hermitian matrix are a phenomenon of codimension 3. This is believed
to be true (although unproved in many cases) for “generic”
Schrödinger operators <cit.>: for a generic
V, the band structure is only expected to show crossings in
codimension 3 (isolated points in 3D, and no crossings in lower
dimensions). The free electron gas, possessing all the symmetries a
Schrödinger operator can have, may thus be considered a worst case
system for disentanglement.




 §.§ The 1D case




  §.§.§ One frozen band, two Wannier functions

The dispersion relation of the 1D free electron gas results in a
crossing between the first and second eigenvalues at half-integers
and between the second and third eigenvalues at integers (see Figure
<ref>). Because of the crossing between the first and second
band, any single Wannier function representing the first band has a
maximal decay of 1/r. We therefore consider the problem of finding
two Wannier functions representing the first band
(N_w = 2, N_f = 1). We do not set an outer window but discretize
the band structure using 2L+1 Fourier modes with L=10 (as we will
see, the Wannier functions we find are compactly supported in Fourier
space, and therefore do not depend on the choice of
L ≥ 2).

For this simple system, initializing the gauge randomly allows our
optimization algorithm to converge robustly to the same Wannier
functions, up to a change of sign and a shift by a lattice vector.
It would therefore seem that the global minimizer of the spread is
unique (up to the invariant degrees of freedom described above) and
real. The SCDM
algorithm with the settings μ = 0 and σ = 2 yields Wannier
functions that are visually indistinguishable from the optimized
Wannier functions, and have a very similar spread (2.51 for the SCDM
algorithm as compared to 2.44 for the optimized Wannier functions).



In Figure <ref> we observe that the first band (which is
frozen) is exactly reproduced. Since the Wannier functions are
localized, the Wannier interpolation is very good, and in particular
has no trouble reproducing the crossing. Furthermore, the optimized
Wannier functions are symmetric (they are real, and the second one is
a translate of the first by half a lattice vector).

However, closer inspection reveals that the optimized Wannier
functions are not exponentially localized, and in fact decay as
1/r^2 for large r (Figure <ref>). The origin of
this slow decay of the Wannier functions is the kink that appears at
ξ = ± 3/2 in their Fourier transform as in Figure <ref>,
which also shows that the Wannier functions are compactly
supported on [-3/2,3/2] in Fourier space. This is because only the
first three bands are occupied by the Wannier functions, as can be
seen via inspection of the U matrix. Since band 3 crosses with band
4 at k=1/2 (corresponding to basis functions K = 1 and K=-2),
the continuity of ψ_nk with respect to k means that
U_31(1/2) = U_32(1/2) = 0. But the first derivative is not zero
at 1/2, which in turns creates the kink in Fourier space at
K=1,k=1/2 and K=-2,k=1/2, that is, ξ = ± 3/2.





The fact that the optimized Wannier functions are only weakly localized is surprising
at first glance, because it was proven in <cit.>
that, for isolated bands, maximally-localized Wannier functions are exponentially localized. What is different in this case? The crucial
point in the analysis of <cit.> is that



the spread is similar to a Dirichlet energy, or a H^1 norm,
in k space. Then the ψ_nk are shown to satisfy an elliptic
equation in k space, which by a bootstrap argument implies their
analyticity. Here, this argument breaks down because the constraint
P_w(k)P_f(k) = P_f(k) is discontinuous at crossings. This
creates an effective “boundary condition” for the gauge U(k) at
crossings that destroys the regularity. A simple analogy is that
eigenfunctions of the Laplacian on [0,1] (critical points of the
Dirichlet energy) are smooth on  when periodic boundary conditions
are imposed, but generically produce kinks at 0 and 1 when
Dirichlet boundary conditions are imposed. Here also, the effective
boundary condition U_31(1/2)=U_32(1/2) = 0 produces a kink.

To remedy this, we show on this particular example how to build
Wannier functions that are of class C^∞ in Fourier space (and therefore
decay faster than any inverse polynomial in real space). 





In order to do so, we define the function α(x) = e^-1/x for
x ≥ 0, α(x) = 0 for x < 0. This function is
C^∞, and identically zero for x ≤ 0. The function
f(x) = α(x)/α(x)+α(1-x) is therefore
C^∞ on ℝ, equal to 0 for x ≤ 0, and equal to
1 for x ≥ 1.

Given U(0^+), U(1/2^±), U(1^-) obtained from the optimized
Wannier functions, we construct a new gauge

    U(k) =
      
        (1-f(2k)) U(0^+) + f(2k) U(1/2^-)     if 0 < k ≤ 1/2,
    
        (1-f(2k-1)) U(1/2^+) + f(2k-1) U(1^-)     if 1/2 < k < 1.

This produces a new set of ψ_nk that are smooth with
respect to k, but not orthogonal and do not span the frozen bands.
We now impose those conditions in the same way we do for the SCDM
procedure as described in Section <ref>. In this specific example, this produces a
smooth gauge, as illustrated in Figure <ref>.



The Wannier functions obtained in this way display more rugged
variation in Fourier space (a consequence of the use of the function
f above), and accordingly the decay is slower for small values of
x. However, because the gauge is smooth, the Wannier functions decay
much faster for large values. Since the gauge is of class C^∞
but not analytic, the Wannier functions decay faster than any
polynomial, but not exponentially. Since the error of Wannier
interpolation on the first band is determined by the interaction of
the Wannier functions on the supercell with their periodic images, the
faster asymptotic decay leads to a better reconstruction of the first
band in this case. Numerical tests indicate that the cross-over point
(above which the reconstruction error on the first band with the
smoothing procedure is smaller than that with the optimized Wannier
functions) is around N=12.



  §.§.§ General case in 1D

In the one-dimensional free electron gas crossings only happen at
half-integers, and between two bands at the same time. Accordingly,
whenever N_w≥ N_f+1 we can construct optimized Wannier functions similar to the
ones above—they are translates of each other, compactly
supported in Fourier space, and decay as 1/r^2. When
N_w = N_f (that is, we are treating a metal as if it was an
insulator), the gauge is discontinuous, and the corresponding Wannier
functions decay as 1/r.



 §.§ The 2D case

In 2D, the first four bands of the free electron gas are degenerate at = (1/2,1/2),
corresponding to the wave vectors 

    = (0,0), (-1,0),(0,-1),(-1,-1).

This means that P_w(1/2,1/2) must span the four-dimensional
subspace corresponding to those four wave vectors. Therefore, freezing
the first band can only produce localized Wannier functions when
N_w≥ 4. Accordingly, we consider the case
N_f = 1 and N_w = 4.

Similarly to the 1D case, the optimized Wannier functions are real,
and differ from each other only by a change of origin. However, they
are not compactly supported in Fourier space, and instead have
decaying components on arbitrarily large wave vectors  (see
Figure <ref>). Their support in Fourier space displays
a checkerboard pattern, and in particular has corners. Furthermore, we
observe numerically that the gradient in Fourier space blows up
near those corners when the size of the k-point grid N is
increased (see Figure <ref>). The maximum value of the
gradient is found to behave as N^1/3. This is consistent with the
singularity of the first eigenfunction of the Laplace operator with
Dirichlet boundary conditions on a domain with corners, which behave
as r^2/3sin(2θ/3) near a corner described by
θ∈ [0,3π/2] in polar coordinates <cit.>.
In particular, their gradient behaves as r^-1/3 near the corner,
explaining the divergence as N^1/3 when discretized on a grid. As
in the 1D case, imposing frozen bands acts as an effective boundary
condition for the gauge and destroys the regularity.

The decay of the optimized Wannier functions in Fourier space is
dictated by the singularity in Fourier space. More precisely, the
first derivative is discontinuous along edges, which corresponds to a
decay of the Wannier functions in the x and y directions as
1/r^2.





 §.§ Discussion

These results show that the maximally-localized Wannier functions only
decay algebraically in general. Although numerical results are harder
to obtain for real materials, we expect our analysis to carry through:
an eigenvalue crossing at a particular  point acts as a constraint
on the gauge, which must at that point be able to span the crossing
eigenspace for the gauge to be continuous. Minimizing the second
moments of the Wannier functions yields a gauge with a
square-integrable but discontinuous first derivative at the crossing
points, resulting in a weak localization.

Our findings are to be contrasted with the recent theoretical result
of <cit.>, which proves under generic hypotheses that there
exists almost-exponentially localized Wannier functions. This simply
means that, unlike in the case of insulators, for entangled band
structures minimizing the second moments is not an asymptotically
optimal strategy.

To get a faster asymptotic decay, one could minimize higher moments.
This corresponds to minimizing integrals of higher derivatives, which
have to be approximated by more complex stencils, and require the
computation of additional overlaps between the u_n than simply
the nearest neighbors. This becomes numerically expensive and complex
to implement. Another possibility is to perform local “smoothing
surgeries” similar to the one we demonstrated in one dimension. This
is likely to be useful only for very fine -point grids however.





§ CONCLUSION AND DISCUSSION


We have developed a variational formulation that, paired with a
specific initialization strategy, is able to robustly construct
Wannier functions for systems with entangled band structure.
Importantly, the definition of Wannier functions must be generalized,
allowing them to lie in a subspace that contains, but is larger than,
the spectral subspace of interest. While this condition adds extra
constraints to our optimization problem, and they can be phrased in
many theoretically equivalent ways, we find one that is particularly
convenient. This results in a formulation that matches that of partly
occupied Wannier functions <cit.>, and
allows us to view the widely used disentanglement procedure as an
alternating minimization algorithm—albeit one that only
takes a single alternation step. As the underlying problem is
non-convex, our choice of initialization strategy via the SCDM
methodology is key. As demonstrated with several real materials, our
method is robust and effective at finding localized functions and
enabling good quality band interpolation. Our variational formulation
is versatile, and can be modified relatively easily to accommodate
additional constraints, such as symmetry constraints, for certain type
of real materials. It would also be interesting to study the behavior
of localization properties of generalized Wannier functions for
systems with non-trivial topological characters.

We also study the free electron gas, providing interesting insights
into the further theoretical study of the localization properties of
generalized Wannier functions. We find that the minimization of the
second moments of the Wannier functions only imposes a relatively weak
algebraic decay. Our analysis suggests that, for real 3D materials,
the disentangled Wannier functions decay asymptotically slowly as
well. Further localization is possible, but the method we present here
is likely to only be useful for very fine -point grids. The
computation of Wannier functions that are localized in both the
pre-asymptotic and asymptotic regime remains an interesting open question.




§ ACKNOWLEDGMENTS


This work was partially supported by the National Science Foundation
under Grant No. DMS-1652330, the Department of Energy under Grants No.
DE-SC0017867 and No. DE-AC02-05CH11231, and the SciDAC project (L. L.);
and the National Science Foundation under grant number DMS-1606277 (A.
D.).






§ GRADIENT OF THE OBJECTIVE FUNCTION





If F is a function from a complex Hilbert space to , we recall that its
gradient (sometimes written ∂ F / ∂ z^*) is defined as the (unique) vector g such that

    F(z+h) = F(z) +  g, h + O(h^2).

All the derivatives and gradient below are in this sense. Note that
this is different from the notion of derivative of a →
function, which is not relevant here (non-trivial complex-to-real
functions are not complex differentiable). The advantage of this
definition is that it allows a straightforward translation of
first-order (but not second-order) optimization algorithms in complex
variables.








We use the numerical setup of <cit.>. Recall
that the -point grid is discretized, with a total of N_
points, and on each  in that grid, we are given a set of
N_g× N_o matrices u_ representing the N_o
orbitals in the outer window discretized on a space of dimension
N_g. The vectors {} are displacements from one -point to a
set of neighbors, and w_ are weights chosen to satisfy

    ∑_ w_^T = I_3,

so that the gradient of a function f() can be approximated
by

    ∇ f() ≈∑_ w_ (f(+)-f())


We look for a set of N_o× N_w matrices {U_} with
orthogonal columns, which define Wannier functions by (<ref>). Let

    M_,,m,n = ⟨ u_m,, u_n,+⟩

be the N_o× N_o overlap matrix between the bands, which is
an input of the algorithm. Then the N_w× N_w overlap matrix between the Wannier
functions defined by U() is

    N_, = U_^*M_,U_+.

The Marzari-Vanderbilt spread functional is given by

    Ω   = ∑_n ||^2_n - |_n|^2,  where
    ^2_n   = 1/N_∑_, w_( 1 - |N_ nn|^2 + (ln N_ nn)^2)
    _n   = - 1/N_∑_, w_ln N_ nn

(equations (11), (31) and (32) of
<cit.>).


We need to compute Ω(U + Δ U) to first
order in Δ U, from which we will identify ∇Ω by
Ω(U+Δ U) - Ω(U) = ∑_((∇Ω)_^* (Δ U)_) + O(Δ U^2). 



We begin with ∑_n⟨ ||^2⟩_n and consider the
following quantity:

    I = ∑_ n f(N_ nn)

where f: →. Then, using the fact that the set of vectors
 is symmetric (contains  as well as -), and that
w_-=w_,

    Δ I    = ∑_ n f”(N_ nn)^*Δ N_ nn
    
         = ∑_ n w_f'(N_ nn)^* (Δ A_^* M_ A_+)_nn + f'(N_ nn)(Δ A_+^*M_^*A_)_nn
       = ∑_ n w_f'(N_ nn)^* (Δ A_^* M_ A_+)_nn + f'(N_-,,nn)(Δ A_^*M_-,^*A_-)_nn
       = ∑_ n w_f'(N_ nn)^* (Δ A_^* M_ A_+)_nn + f'(N_,-,nn^*)(Δ A_^*M_,-^*A_-)_nn
       = ∑_ n w_(f'(N_ nn)^* + f'(N_ nn^*)) (Δ A_^* M_ A_+)_nn

and the gradient is therefore

    (∇ I)_ mn   = ∑_ w_(f'(N_ nn)^* + f'(N_ nn^*)) (M_ A_+)_mn

It can be checked using similar arguments that, when

    I =∑_ n
    g(N_ nn),

then

    (∇ I)_ mn   = -∑_ w_(g'(N_ nn)^* + g'(N_ nn^*)) (M_ A_+)_mn

Applying these formulas with

    f(z)    = 1-|z|^2 + (ln z)^2,   f'(z) = -2z + 2 iln z/z^*
    
      g(z)    = -ln z,    g'(z) = -i/z^*

we get

    (∇Ω)_ mn   = 4/N_∑_(-N_ nn^* - iln N_ nn +  r
      _n·/N_ nn)(M_ A_+)_mn


Note that this is the unconstrained gradient of Ω with
respect to U. Using the chain rule one with U()=[ I_N_f     0;     0     Y ] X one can easily
derive the gradient of Ω with respect to (X,Y). Then we
simply have to minimize with respect to (X,Y) subject to the
orthogonality constraints for X and Y using standard methods <cit.>.
















siam


