---
abstract: 'We consider a variational problem for three-dimensional (3D) classical lattice models. We construct the trial state as a two-dimensional product of local variational weights that contain auxiliary variables. We propose a stable numerical algorithm for the maximization of the variational partition function per layer. The numerical stability and efficiency of the new method are examined through its application to the 3D Ising model.'
author:
- |
    Andrej <span style="font-variant:small-caps;">Gendiar</span>,$^{1,2}$ Nobuya <span style="font-variant:small-caps;">Maeshima</span>,$^3$ and Tomotoshi <span style="font-variant:small-caps;">Nishino</span>$^1$\
    $^1$Department of Physics, Faculty of Science, Kobe University,\
    Kobe 657-8501, Japan\
    $^2$Institute of Electrical Engineering, Slovak Academy of Sciences,\
    Dúbravská cesta 9, SK-842 39 Bratislava, Slovakia\
    $^3$Department of Physics, Graduate School of Science, Osaka University,\
    Toyonaka 560-0043, Japan
title: Stable Optimization of a Tensor Product Variational State
---

Introduction
============

The density matrix renormalization group (DMRG) has been applied to a variety of one-dimensional (1D) quantum systems and 2D classical systems [@White; @Nishino; @DMRG]. This method has also been applied to finite-size 2D quantum systems, which can be represented as 1D quantum systems with long-range interactions. [@Liang; @DMRG] Despite this success, no extension of the DMRG to [*infinitely large*]{} higher-dimensional systems has been reported to this time. This is partially because the density matrix eigenvalues decay very slowly in higher-dimensional systems, and the RG transformation, [*which is obtained from the diagonalization of the density matrix*]{}, becomes ineffective. [@Peschel]

The numerical efficiency of the DMRG for 1D quantum and 2D classical systems is due in part to its variational background, [@Ost; @Rom; @Sierra] in which the trial state is constructed as a product of orthogonal matrices that represent the block-spin (or the renormalization group) transformations. Such a construction of the variational (or trial) state can be generalized to higher dimensions. A simple example is to use a 2D classical system as a variational state for 2D quantum and 3D classical systems. Martín-Delgado [*et al.*]{} employed the 6-vertex model as a trial wave function for 2D lattice spin/electron systems. [@Sierra_st] Okunishi and Nishino considered the direct extension of the Kramers-Wannier approximation [@KWA] to the 3D Ising model, in which the variational state is the 2D Ising model subject to an effective magnetic field. [@KWON] These examples have demonstrated the potential usefulness of the 2D variational state, constructed as a product of local variational weights.

For the purpose of obtaining a better variational result, Nishino [*et al.*]{} developed a numerical method, the tensor product variational approach (TPVA), which automatically optimizes the 2D variational state using the solution of a self-consistent equation. [@TPVA1] They applied this method to both the 3-state ($q = 3$) Potts and the Ising models on simple cubic lattices, and reported that TPVA yields a better estimate of the transition temperature for the $q = 3$ Potts model. This is because the number of variational parameters in the trial state for the Potts model, $3^4_{~}$, is greater than that for the Ising model, $2^4_{~}$.

In general, the TPVA gives a lower variational free energy when the trial state contains a larger number of parameters. A way of increasing this number is to employ a variational state that contains auxiliary variables, for example block spin variables. [@TPVA2][^1] This generalization, however, introduces a serious instability into the numerical optimization of the variational state. [@TPVA2]

In this paper we report on a method of stabilizing the numerical optimization process. We introduce an orthogonal condition placed on a small change of the local variational weight. This condition prevents an ‘unexpected decrease’ of the norm of the variational state, which caused numerical instability in previous calculations. [@TPVA2] In the next section, we review the formulation of the TPVA. [@TPVA2; @TPVA1; @TPVA3] In §3 we consider the stability of the optimization process. The numerical efficiency of the stabilized TPVA is examined in §4 in the case that it is applied to the 3D Ising model. Conclusions are given in §5.

Tensor product variational approach
===================================

![The Boltzmann weight $W_{\rm B}^{(ij)}\{ \bar\sigma | \sigma \}$ and the transfer matrix ${\cal T}[ \bar\sigma | \sigma ]$ in the case $N=2$.[]{data-label="tm"}](tm.eps){width="85mm"}

We consider the 3D Ising model on a simple cubic lattice as an example. The system size is $2N \times 2N \times \infty$ in the $X$-, $Y$-, and $Z$-directions. On each lattice point there is an Ising spin $\sigma$, with $\sigma = \pm 1$. A ferromagnetic interaction $-J \sigma \sigma'$ exists between nearest-neighbor spins.

As shown in Fig. \[tm\], the transfer matrix ${\cal T}$ connects adjacent spin layers $[ \sigma ]$ and $[ \bar\sigma ]$, where $[ \sigma ]$ represents all the spins in a layer: $$[ \sigma ] = \{ \sigma_{1\,1}^{~}, \sigma_{2\,1}^{~}, \cdots,
\sigma_{2N\,1}^{~}, \sigma_{1\,2}^{~}, \sigma_{2\,2}^{~}, \cdots,
\sigma_{2N\,2N}^{~} \} \,.
\label{sm}$$ For simplicity, we consider a symmetric transfer matrix, $${\cal T}[ \bar\sigma | \sigma ] = \prod_{i,j = 1}^{2N - 1}
W_{\rm B}^{(ij)}\{ \bar\sigma | \sigma \} \, ,
\label{irfeqtm}$$ where $W_{\rm B}^{(ij)}\{ \bar\sigma | \sigma \}$ represents a local Boltzmann weight with respect to a unit cube at the position $(X,Y)=(i,j)$ (see Fig. \[tm\]). We denote a spin plaquett, a group of the nearest 4 spins, as $$\{ \sigma \}
= (\sigma_{i j}^{~}\ \sigma_{i j+1}^{~}
\ \sigma_{i+1 j}^{~}\ \sigma_{i+1 j+1}^{~})
= (\sigma_{ij}^{~}\ \sigma_{ij^\prime}^{~}
\ \sigma_{i^\prime j}^{~}\ \sigma_{i^\prime j^\prime}^{~}),$$ using the $i^\prime = i + 1$ and $j^\prime = j + 1$. With this notation, we can write the local Boltzmann weight of the 3D Ising model as follows: $$\begin{aligned}
\nonumber
W_{\rm B}^{(ij)}\{ \bar\sigma | \sigma \}
& = & \exp \biggl[-\frac{J}{4 k_{\rm B}^{~} T} \left(
\sigma_{ij}^{~}\sigma_{i'j}^{~}+\sigma_{i'j}^{~}\sigma_{i'j'}^{~}
+\sigma_{i'j'}^{~}\sigma_{ij'}^{~}
+\sigma_{ij'}^{~}\sigma_{ij}^{~}+\bar\sigma_{ij}^{~}\bar\sigma_{i'j}^{~}
\right. \\
\label{WB3D}
\nonumber
& + &
\bar\sigma_{i'j}^{~}\bar\sigma_{i'j'}^{~}
+\bar\sigma_{i'j'}^{~}\bar\sigma_{ij'}^{~}+
\bar\sigma_{ij'}^{~}\bar\sigma_{ij}^{~}\left.
+\sigma_{ij}^{~}\bar\sigma_{ij}^{~}+
\sigma_{i'j}^{~}\bar\sigma_{i'j}^{~}+
\sigma_{i'j'}^{~}\bar\sigma_{i'j'}^{~}+
\sigma_{ij'}^{~}\bar\sigma_{ij'}^{~}
\right)\biggr] \, .\end{aligned}$$ We have thus expressed the 3D Ising model as a special case of the ‘interaction round a face’ (IRF) model.

![A graphical representation of the local tensor $V^{(ij)}_{~}$ and the construction of the trial state $\Psi$ in the case $N = 2$. The circles and ovals denote the Ising spins $\sigma$ and the $n$-state auxiliary variables $\xi$, respectively.[]{data-label="vs"}](vs.eps){width="85mm"}

For an arbitrary variational function $\Psi$, the Rayleigh ratio $$\lambda( \Psi ) = \frac{\sum\limits_{[ \bar\sigma ], [ \sigma ]}
\Psi\left[ \bar\sigma \right] {\cal T}[ \bar\sigma | \sigma ]
\Psi\left[ \sigma\right ]}{\sum\limits_{[ \bar\sigma ], [ \sigma ]}
\Psi\left[ \bar\sigma\right ] \Psi\left[ \sigma\right ]}
\equiv\frac{\langle \Psi | {\cal T} | \Psi \rangle}{\langle \Psi |
\Psi \rangle}
\label{rr}$$ gives the variational partition function per layer. In the framework of the TPVA, the trial state $\Psi$ is constructed as a contracted product of local variational weights as $$\Psi\left[ \sigma\right ]
\, = \, \sum\limits_{[\xi]}^{~}\prod\limits_{i,j=1}^{2N-1}
V\left( {{\sigma_{ij}^{~}}\atop{\xi_{ij}^{~}}}\ {{\sigma_{i^\prime
j}^{~}}\atop{\xi_{i^\prime j}^{~}}}
\ {{\sigma_{ij^\prime}^{~}}\atop{\xi_{ij^\prime}^{~}}}\
{{\sigma_{i^\prime j^\prime}^{~}}
\atop{\xi_{i^\prime j^\prime}^{~}}} \right)
\, \equiv \, \sum\limits_{[ \xi ]}
\prod\limits_{i,j=1}^{2N-1}V^{(ij)}\left( {{\{ \sigma \}}\atop{\{ \xi
\}}} \right) \, ,
\label{vsIRF}$$ where the variables $\xi_{ij}^{~}$ of the local weight $V$ are auxiliary variables that can be in one of $n$ states ($0,1,\ldots,n-1$), while the spin variables $\sigma_{ij}^{~}$ can be in one of 2 states ($\pm 1$) [^2] (see Fig. \[vs\]). We have expressed a group of auxiliary variables as $[ \xi ]$, analogously to $[ \sigma ]$ in Eq. (\[sm\]). [^3] We are interested in the bulk properties of lattice models, and therefore we consider the case in which the system size $2N$ is sufficiently large and assume that the local variational weight $V^{(ij)}_{~}$ is position independent. Hereafter, we refer to this variational state as the ‘tensor product state’ (TPS).

Because both the trial state $\Psi$ and the transfer matrix ${\cal T}$ are written as products of local variational weights, both the numerator and denominator of Eq. (\[rr\]) can also be expressed as products of stacked local weights: $$\begin{aligned}
\langle \Psi | {\cal T} | \Psi \rangle & = & \sum\limits_{{[
\bar\sigma ],[ \sigma ]}
\atop{[ \bar\xi ], [ \xi ]}}\prod\limits_{i,j=1}^{2N-1}
V^{(ij)}_{~}\left( {{\{ \bar\sigma \}}\atop{\{ \bar\xi \}}} \right)
W_{\rm B}^{(ij)}\{ \bar\sigma | \sigma \}\
V^{(ij)}_{~}\left( {{\{ \sigma \}}\atop{\{ \xi \}}} \right) \, , \nonumber\\
\langle \Psi | \Psi \rangle & = & \sum\limits_{{[ \sigma ]}
\atop{[ \bar\xi ], [ \xi ]}}\prod\limits_{i,j=1}^{2N-1}
V^{(ij)}_{~}\left( {{\{ \sigma \}}\atop{\{ \bar\xi \}}} \right)\
V^{(ij)}_{~}\left( {{\{ \sigma \}}\atop{\{ \xi \}}} \right) \, .
\label{vv}\end{aligned}$$ In other words, both $\langle \Psi | \Psi \rangle$ and $\langle \Psi |
{\cal T} | \Psi \rangle$ are partition functions of stacked two-dimensional classical systems. This structure enables us to calculate $\langle \Psi | \Psi \rangle$ and $\langle \Psi | {\cal T} | \Psi \rangle$ accurately by use of numerical renormalization techniques. [@CTMRG1; @CTMRG2; @Nishino] Thus, the variational partition function $\lambda( \Psi )$ can be easily calculated for an arbitrary TPS constructed from the local variational weight $V$.

![Graphical representation of the ‘matrices’ ${\cal A}$ (on the left) and ${\cal B}$ (on the right) in Eqs. (\[eqa\]) and (\[eqb\]), respectively.[]{data-label="ab"}](ab.eps){width="130mm"}

We maximize $\lambda( \Psi )$ by tuning elements of the local weight $V$. In order to clarify the variational structure with respect to $V$, let us divide $\langle \Psi | \Psi \rangle$ into two parts: ([**i**]{}) the adjacent local weights $V^{(NN)}_{~}$ and $\bar V^{(NN)}_{~}$ at the center, and ([**ii**]{}) the rest, consisting of the stacked classical system with a puncture at the center [@Martin], expressed as $${\cal A}\left( {{\{ \sigma \}}\atop{\{ \bar\xi \}}}\left | {{\{
\sigma \}}\atop{\{ \xi \}}}
\right. \right) = \sum\limits_{{[ \sigma ]^\prime}\atop{[ \bar\xi
]^\prime, [ \xi ]^\prime}}
\prod\limits_{(i,j)\neq(N,N)}V^{(ij)}_{~}\left( {{\{ \sigma \}}
\atop{\{ \bar\xi \}}} \right)V^{(ij)}_{~}\left( {{\{ \sigma \}}\atop{\{
\xi \}}} \right) \, ,
\label{eqa}$$ where the configuration sums are taken over all values of the variables $\sigma$ and $\xi$, except for those at the center that are variables of $V^{(NN)}_{~}$ and $\bar V^{(NN)}_{~}$. The notation $\prod_{(i,j)\neq(N,N)}$ indicates that $V^{(NN)}_{~}$ is not included in the product. In the same manner, we divide $\langle \Psi | {\cal T} |\Psi \rangle$ into ([**i’**]{}) $V^{(NN)}_{~}$ and $\bar V^{(NN)}_{~}$, and ([**iii**]{}) the rest, expressed by $$\begin{aligned}
\nonumber
{\cal B}\left({{\{ \bar\sigma \}}\atop{\{ \bar\xi \}}}\left |
{{\{ \sigma \}}\atop{\{ \xi \}}}\right.\right)
& = &
W_{\rm B}^{(NN)}\{ \bar\sigma | \sigma \}\sum\limits_{{[ \bar\sigma ]^\prime,
[ \sigma ]^\prime}\atop{[ \bar\xi ]^\prime, [ \xi
]^\prime}}\prod\limits_{(i,j)\neq(N,N)}
V^{(ij)}_{~}\left({{\{ \bar\sigma \}}\atop{\{ \bar\xi \}}}\right)\\
&\times&W_{\rm B}^{(ij)}\{ \bar\sigma | \sigma \}V^{(ij)}_{~}\left({{\{ \sigma \}}
\atop{\{ \xi \}}}\right),
\label{eqb}\end{aligned}$$ which corresponds to a partially punctured stacked classical system. The meanings of the matrices ${\cal A}$ and ${\cal B}$ are represented graphically in Fig. \[ab\]. Using this new notation, we can rewrite the Rayleigh ratio, Eq. (\[rr\]), as the following ratio of bilinear forms: $$\label{vbvvav}
\lambda( \Psi )
  \, = \,
\frac{\sum\limits_{{\{ \bar\sigma \}, \{ \sigma \}}\atop{\{ \bar\xi
\},\{ \xi \}}}
V^{(NN)}_{~}\left({{\{ \bar\sigma \}}\atop{\{ \bar\xi \}}}\right)
{\cal B}\left({{\{ \bar\sigma \}}\atop{\{ \bar\xi \}}}
\left | {{\{ \sigma \}}\atop{\{ \xi
\}}}\right.\right)V^{(NN)}_{~}\left({{\{ \sigma \}}
\atop{\{ \xi \}}}\right)}
{\sum\limits_{{\{ \sigma \}}\atop{\{ \bar\xi \},\{ \xi \}}}
V^{(NN)}_{~}\left({{\{ \sigma \}}\atop{\{ \bar\xi \}}}\right)
{\cal A}\left({{\{ \sigma \}}\atop{\{ \bar\xi \}}}
\left|{{\{ \sigma \}}\atop{\{ \xi
\}}}\right.\right)V^{(NN)}_{~}\left({{\{ \sigma \}}
\atop{\{ \xi \}}}\right)}
\, \equiv \, \frac{( V | {\cal B} | V )}{( V | {\cal A} | V )} \, .$$ Here, we have interpreted the variational weight at the center $V^{(NN)}_{~}$ as a $(2n)^4_{~}$-dimensional vector and have written it simply as $| V )$.

As we have assumed that the system size $2N$ is sufficiently large, the variation of $\lambda( \Psi )$ with respect to a uniform modification of local weights, $$\delta\lambda( \Psi ) = \sum_{i j}^{~}
\frac{\partial \lambda( \psi )}{\partial V^{(ij)}_{~}}
\delta V^{(ij)}_{~}+ {\rm [higher~order~terms]},$$ is almost proportional to $\partial \lambda( \psi )
/ \partial V^{(NN)}_{~}$, which is the contribution from the local variation only at the center. [@TPVA1; @TPVA2] From the optimal condition $\partial\lambda( \psi ) / \partial V^{(NN)}_{~} = 0$, we obtain the generalized eigenvalue problem $${\cal B} \, | V ) \, = \, \lambda {\cal A} \, | V ) \, .
\label{sce}$$ We use this equation when we optimize the TPS. Note that Eq. (\[sce\]) is a non-linear equation with respect to the local weights $V$, since the $(2n)^4_{~}$-dimensional matrices ${\cal A}$ and ${\cal B}$ themselves are functionals of $V^{(NN)}_{~}$. Therefore, Eq. (\[sce\]) should be solved self-consistently by use of successive, small improvements of the local variational weight $V^{(NN)}_{~}$.

In Ref. [@TPVA2], Nishino [*et al.*]{} apply the TPVA to the 3D Ising model, which is represented as a symmetric 16-vertex model. For the case $n = 2$, they succeeded in optimizing the TPS. However, they encountered numerical instability in the case $n \geq 3$. The reason for this is that the matrix ${\cal A}$ is not always positive definite during the process in which the variational state is being improved, although ${\cal A}$ should be positive definite after the optimization is completed. We discuss the cause of this instability and propose a method of stabilization in the following.

Stabilization of the self-consistent improvement process
========================================================

Consider an infinitesimal change of the local weight at the center of the system, $$| V ) \, \rightarrow \, | V' ) \, = \, | V ) \, + \, \varepsilon | X ),
\label{VX}$$ where $| X )$ is an arbitrary $(2n)^4_{~}$-dimensional vector. The Rayleigh ratio, Eq. (\[vbvvav\]), is modified as $$\lambda' \, \equiv \, \frac{( V' | {\cal B} | V' )}{( V' | {\cal A} | V' )}
               \, = \, \frac{( V  | {\cal B} | V  )}{( V  | {\cal A} | V  )}
\left(
\frac{\displaystyle 1 + 2 \varepsilon\frac{( V | {\cal B} | X )}{( V |
{\cal B} | V )}
+ \varepsilon^2\frac{( X | {\cal B} | X )}{( V | {\cal B} | V )}}
{\displaystyle 1 + 2\varepsilon \frac{( V |{\cal A} | X )}{( V |
{\cal A} | V )}
+ \varepsilon^2\frac{( X | {\cal A} | X )}{( V | {\cal A} | V )}}\right) \, .
\label{der1}$$ If $\varepsilon$ is sufficiently small, we can expand $\lambda'$ as $$\lambda' = \lambda \left[ 1 + 2 \varepsilon ( \beta - \alpha )
  + {\cal O}(\varepsilon^2)\right],
\label{der2}$$ with $$\lambda = \frac{( V | {\cal B} | V )}{( V | {\cal A} | V )} \, ,\quad
\alpha  = \frac{( V | {\cal A} | X )}{( V | {\cal A} | V )} \, ,\quad
\beta   = \frac{( V | {\cal B} | X )}{( V | {\cal B} | V )} \, .$$ It seems appropriate to select $| X )$ to maximize $\beta - \alpha$. However, such a choice of $| X )$ tends to cause the expectation value of the denominator ${( V' | {\cal A} | V' )}$ to decrease, and after several self-consistent iterations, ${( V | {\cal A} | V )}$ becomes very small or negative. This is a cause of numerical instability in the previous TPVA algorithm. We prevent this [*anomalous decrease*]{} of $( V | {\cal A} | V )$ by choosing $| X )$ to satisfy $$( V' | {\cal A} | V' ) \, = \, ( V | {\cal A} | V )
                           + {\cal O}(\varepsilon^2_{~}) \, .
\label{shrink}$$ This is equivalent to choosing $| X )$ to be orthogonal to ${\cal A} | V )$, which gives $\alpha = 0$. With such a choice, the maximization of $\lambda$ is carried out with respect to $\beta$ only. A reasonable choice of $| X )$ can be obtained using the Schmidt orthogonalization: $$| X ) = {\cal B} | V ) - {\cal A} | V ) \, ( V | {\cal AB} | V ) \, .
\label{orthog}$$ The $| X )$ thus obtained at least yields non-negative $\beta$. We include $| X )$ in the self-consistent calculation of the TPVA as follows:

-   Prepare an arbitrary initial variational weight $| V )$.

-   Calculate the matrices ${\cal A}$ and ${\cal B}$ using a numerical RG method.

-   Obtain $| X )$ from Eq. (\[orthog\]).

-   Improve the local weight according to $| V^{\prime} ) = | V ) +
    \varepsilon | X )$.

-   Terminate the modification of $| V )$ if the computed thermodynamic functions have converged; otherwise go back to step (2).

The third step is the main difference between the present and the previous TPVA algorithms. Since $( V | {\cal A} | V )$ changes only by an amount of order $\varepsilon^2_{~}$ at each iteration, a large number of iterations is necessary to realize convergence, especially when $\varepsilon$ is small.

Numerical results
=================

We now confirm the numerical stability of the new algorithm in the case that it is applied to the 3D Ising model. Hereafter, we set $J/k_{\rm B}^{~}=-1$ and denote by $m$ the number of the states of the block-spin variable used in the numerical RG calculation [@CTMRG1; @CTMRG2], which is used for the calculation of ${\cal A}$ and ${\cal B}$.

Figure \[Convg\] shows the convergence of the spontaneous magnetization $\langle \sigma \rangle$ with respect to the number of numerical iterations at the temperature $T = 4.504$. (The selected temperature $T=4.504$ is slightly lower than the critical temperature $T_{\rm c}^{~} = 4.5115$ obtained by Monte Carlo simulations [@Mc] which is considered to be the most reliable result.) The parameter $\varepsilon$ in Eq. (\[VX\]) is chosen to be $10^{-3}_{~}$ for normalized $| V )$ and $| X )$. In both the cases $n = 1$ and $n = 2$, $\langle \sigma \rangle$ begins to converge monotonically to the final value after several hundred iterations.

![Convergence of the spontaneous magnetization $\langle \sigma \rangle$ with respect to the number of numerical steps at $T = 4.504$.[]{data-label="Convg"}](Mg_cnvg.eps){width="100mm"}

Figure \[IsingIRF\](a) displays $\langle\sigma\rangle$ thus calculated for $n = 1$ (no auxiliary variables) and $n = 2$ with $m = 5$. In the region $T < 4.2$, the difference between these cases is not visible. As shown in the inset, near the critical temperature, the calculated $\langle \sigma \rangle$ with $n = 2$ decays more rapidly than that with $n = 1$. The estimated transition temperatures, where the obtained $\langle \sigma \rangle$ falls to zero, are $T = 4.57$ ($n = 1$) and 4.55 ($n = 2$).

To this point, we have expressed the 3D Ising model as a 3D IRF model. The Ising model can also be expressed as a special case of the symmetric 64-vertex model. [@TPVA2] Applying the stabilized TPVA algorithm to this vertex expression, we obtain the $\langle \sigma \rangle$ plotted in Fig. \[IsingIRF\](b). In this case, the estimated transition temperatures are $T = 4.533$ (with $n = 2$ and $m = 16$) and 4.525 (with $n = 3$ and $m = 12$). All of these calculated transition temperatures are higher than $T_{\rm c}^{\rm MC} = 4.5115$, obtained by Monte Carlo simulation. [@Mc] It is thus seen that with the TPVA, the ordered phase tends to be stabilized.

![Calculated spontaneous magnetizations. (a) $\langle\sigma\rangle$ obtained for $( n, m ) = ( 2, 5 )$. The curve here was drawn using the Spline interpolation. The inset displays $\langle\sigma\rangle$ near the criticality for $( n, m ) = ( 1, 5 )$ (circles) and $( 2, 5 )$ (triangles). (b) $\langle\sigma\rangle$ obtained from the TPVA applied to the Ising model represented as the 64-vertex model [@TPVA2] for $( n, m ) = ( 2, 16 )$. The inset plots the data for $( n, m ) = ( 2, 16 )$ (circles) and $( n, m ) = ( 3, 12 )$ (triangles).[]{data-label="IsingIRF"}](Mg_IRF.eps "fig:"){width="67mm"} ![Calculated spontaneous magnetizations. (a) $\langle\sigma\rangle$ obtained for $( n, m ) = ( 2, 5 )$. The curve here was drawn using the Spline interpolation. The inset displays $\langle\sigma\rangle$ near the criticality for $( n, m ) = ( 1, 5 )$ (circles) and $( 2, 5 )$ (triangles). (b) $\langle\sigma\rangle$ obtained from the TPVA applied to the Ising model represented as the 64-vertex model [@TPVA2] for $( n, m ) = ( 2, 16 )$. The inset plots the data for $( n, m ) = ( 2, 16 )$ (circles) and $( n, m ) = ( 3, 12 )$ (triangles).[]{data-label="IsingIRF"}](Mg_Vertex.eps "fig:"){width="67mm"}

We finally compare the stabilized numerical algorithm considered here with the previous one. [@TPVA2] The stabilization provided by the Schmidt orthogonalization, expressed by Eq. (\[orthog\]), enables us to perform the calculation for those cases $n = 2$ (IRF expression) and $n = 3$ (Vertex expression) for which the previous algorithm does not yield a convergent result. With regard to the computational time required to obtain the convergent numerical result for the cases $n = 1$ (IRF expression) and $n = 2$ (Vertex expression), for which the previous and the stabilized algorithms give the same numerical result, the stabilized algorithm is about 10 times slower than the previous algorithm. In the variational calculation employing the TPVA, the numerator $( V | {\cal B} | V )$ should be maximized and the denominator $( V | {\cal A} | V )$ should be minimized through successive improvements of $V$. Though $( V | {\cal B} | V )$ increases rather rapidly, $( V | {\cal A} | V )$ decreases slowly, because the stabilization condition Eq. (\[shrink\]) restricts the size of the change of $( V | {\cal A} | V )$ to order $\varepsilon^2_{~}$.

Conclusion and discussion
=========================

We have proposed a stabilized partition function maximization process for numerical calculations employing the TPVA, imposing the orthogonality condition on the small change made to the local variational weight. This improvement enables us to obtain the magnetization using a TPS that contains $n$-state auxiliary variables, in particular, with $n = 2$ for an IRF-type TPS and $n = 3$ for a vertex-type TPS.

The orthogonality condition expressed by Eq. (\[orthog\]) stabilizes the numerical calculation, but it also slows the convergence to the variational maximum, because it causes the improvement of the denominator of the variational formulation to be very slow. Contrastingly, the vertical density matrix approach (VDMA) [@TPS], which also uses a TPS as its trial state, exhibits rapid convergence. In the VDMA, the local weight is improved through a RG transformation obtained from the diagonalization of the density matrix, and thus both the denominator and the numerator of the variational ratio are improved equally rapidly. The VDMA, however, has the shortcoming that it does not fully improve the TPS, because the RG transformation used in the VDMA is not determined so as to maximize the variational partition function per a layer. [*In higher dimensions, the direct diagonalization of the density matrix does not give the most appropriate RG transformation.*]{} Our next project is to combine the advantages of the TPVA and the VDMA.

Acknowledgements {#acknowledgements .unnumbered}
================

We thank K. Okunishi, Y. Hieida and Y. Akutsu for valuable discussions. This work has been partially supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture (Grant No. 09640462 and No. 11640376) and by the Slovak Grant Agencies, VEGA No. 2/7201/21 and 2/3118/23. A. G. is also supported by the Japan Society for the Promotion of Science (P01192).

[99]{} S. R. White, Phys. Rev. Lett. [**69**]{},2863 (1992), Phys. Rev. B [**48**]{}, 10345 (1993). *Density-Matrix Renormalization – A New Numerical Method in Phys.*, Lecture Notes in Physics [**528**]{}, ed. I. Peschel, X. Wang, M. Kaulke and K. Hallberg (Springer Verlag, Berlin, 1999). T. Nishino, J. Phys. Soc. Jpn. [**64**]{}, 3598 (1995). S. D. Liang and H. B. Pang, Phys. Rev. B [**49**]{}, 9214 (1994). M. C. Chung and I. Peschel, Phys. Rev. B [**64**]{}, 064412 (2001). S. Östlund and S. Rommer, Phys. Rev. Lett. [**75**]{}, 3537 (1995). S. Rommer and S. Östlund, Phys. Rev. B [**55**]{}, 2164 (1997). J. Dukelsky, M.A. Martín-Delgado, T. Nishino and G. Sierra, Europhys. Lett. [**43**]{}, 457 (1998). M.A. Martín-Delgado, M. Roncaglia, and G. Sierra, Phys. Rev. B [**64**]{}, 075117 (2001). H.A. Kramers and G.H. Wannier, Phys. Rev. [**60**]{}, 263 (1941). K. Okunishi and T. Nishino, Prog. Theor. Phys. [**103**]{}, 541 (2000). T. Nishino, K. Okunishi, Y. Hieida, N. Maeshima and Y. Akutsu, Nucl. Phys. B [**575**]{}, 504 (2000). T. Nishino, K. Okunishi, Y. Hieida, N. Maeshima, Y. Akutsu, and A. Gendiar, Prog. Theor. Phys. [**105**]{}, 409 (2001). N. Maeshima, Y. Hieida, Y. Akutsu, T. Nishino, and K. Okunishi, Phys. Rev. E [**64**]{}, 106705 (2001). A. Gendiar and T. Nishino, Phys. Rev. E [**65**]{}, 046702 (2002). T. Nishino, K. Okunishi, J. Phys. Soc. Jpn. [**65**]{}, 891 (1996). T. Nishino, K. Okunishi, J. Phys. Soc. Jpn. [**66**]{}, 3040 (1997). M.A. Martín-Delgado, J. Rodriguez-Laguna, and G. Sierra, Nucl. Phys. B [**601**]{}, 569 (2001). W. Janke and R. Villanova, Nucl. Phys. B [**489**]{}, 679 (1997).

[^1]: There is another way of tuning the local variational weight by way of the density matrix renormalization in vertical direction (see Ref. [@TPS]).

[^2]: Two kinds of tensor products are known. One is the IRF type [@TPVA1] and the other is the vertex type. [@TPVA2] We use the former throughout this section.

[^3]: When $n = 1$, the variational function $\Psi$ does not contain any auxiliary variables.
