# DESY 99-170 Least-squares optimized polynomials for fermion simulations11footnote 1Talk given at the workshop on Numerical Challenges in Lattice QCD, August 1999, Wuppertal University; to appear in the proceedings.

## 1 Introduction

In popular Monte Carlo simulation algorithms for QCD and other similar quantum field theories the main difficulty is the evaluation of the determinant of the fermion action matrix. This can be achieved by stochastic procedures with the help of auxiliary bosonic “pseudofermion” fields.

In the two-step multi-bosonic (TSMB) algorithm  an approximation of the fermion determinant is achieved by the pseudofermion fields corresponding to a polynomial approximation of some negative power $`x^\alpha `$ of the fermion matrix . The auxiliary bosonic fields are updated according to the multi-bosonic action . The error of the polynomial approximation is corrected in a global accept-reject decision by using better polynomial approximations. Sometimes a reweighting of gauge configurations in the evaluation of expectation values is also necessary. This can also be performed by high order polynomials.

The polynomials used in the TSMB algorithm have to approximate the function $`x^\alpha \overline{P}(x)`$ in some non-negative interval $`x[ϵ,\lambda ],\mathrm{\hspace{0.33em}0}ϵ<\lambda `$. Here $`\overline{P}(x)`$ is a known polynomial, typically another cruder approximation of $`x^\alpha `$. The approximation scheme and optimization procedure can be chosen differently. The least-squares optimization  is an efficient and flexible possibility. (For other approximation schemes see .)

In this review the basic relations for least-squares optimized polynomials are presented as introduced in . Particular attention is paid to a recurrence scheme which can be applied for determining the necessary high order polynomials and for evaluating them numerically. The details of the TSMB algorithm will not be considered. For a comprehensive summary and references see . For experience on the application of TSMB in a recent large scale numerical simulation see .

## 2 Basic relations

Least-squares optimization provides a general and flexible framework for obtaining the necessary optimized polynomials in multi-bosonic fermion algorithms. Here we introduce the basic formulae in the way it has been done in .

We want to approximate the real function $`f(x)`$ in the interval $`x[ϵ,\lambda ]`$ by a polynomial $`P_n(x)`$ of degree $`n`$. The aim is to minimize the deviation norm

$$\delta _n\left\{N_{ϵ,\lambda }^1_ϵ^\lambda 𝑑xw(x)^2\left[f(x)P_n(x)\right]^2\right\}^{\frac{1}{2}}.$$
(1)   
Here $`w(x)`$ is an arbitrary real weight function and the overall normalization factor $`N_{ϵ,\lambda }`$ can be chosen by convenience, for instance, as

$$N_{ϵ,\lambda }_ϵ^\lambda 𝑑xw(x)^2f(x)^2.$$
(2)   
A typical example of functions to be approximated is $`f(x)=x^\alpha /\overline{P}(x)`$ with $`\alpha >0`$ and some polynomial $`\overline{P}(x)`$. The interval is usually such that $`0ϵ<\lambda `$. For optimizing the relative deviation one takes a weight function $`w(x)=f(x)^1`$.

$`\delta _n^2`$ is a quadratic form in the coefficients of the polynomial which can be straightforwardly minimized. Let us now consider, for simplicity, only the relative deviation from the simple function $`f(x)=x^\alpha =w(x)^1`$. Let us denote the polynomial corresponding to the minimum of $`\delta _n`$ by

$$P_n(\alpha ;ϵ,\lambda ;x)\underset{\nu =0}{\overset{n}{}}c_{n\nu }(\alpha ;ϵ,\lambda )x^{n\nu }.$$
(3)   
Performing the integral in $`\delta _n^2`$ term by term we obtain

$$\delta _n^2=12\underset{\nu =0}{\overset{n}{}}c_\nu V_\nu ^{(\alpha )}+\underset{\nu _1,\nu _2=0}{\overset{n}{}}c_{\nu _1}M_{\nu _1,\nu _2}^{(\alpha )}c_{\nu _2},$$
(4)   
where

$$V_\nu ^{(\alpha )}=\frac{\lambda ^{1+\alpha +n\nu }ϵ^{1+\alpha +n\nu }}{(\lambda ϵ)(1+\alpha +n\nu )},$$

$$M_{\nu _1,\nu _2}^{(\alpha )}=\frac{\lambda ^{1+2\alpha +2n\nu _1\nu _2}ϵ^{1+2\alpha +2n\nu _1\nu _2}}{(\lambda ϵ)(1+2\alpha +2n\nu _1\nu _2)}.$$
(5)   
The coefficients of the polynomial corresponding to the minimum of $`\delta _n^2`$, or of $`\delta _n`$, are

$$c_\nu c_{n\nu }(\alpha ;ϵ,\lambda )=\underset{\nu _1=0}{\overset{n}{}}M_{\nu \nu _1}^{(\alpha )1}V_{\nu _1}^{(\alpha )}.$$
(6)   
The value at the minimum is

$$\delta _n^2\delta _n^2(\alpha ;ϵ,\lambda )=1\underset{\nu _1,\nu _2=0}{\overset{n}{}}V_{\nu _1}^{(\alpha )}M_{\nu _1,\nu _2}^{(\alpha )1}V_{\nu _2}^{(\alpha )}.$$
(7)   
The solution of the quadratic optimization in (6)-(7) gives in principle a simple way to find the required least-squares optimized polynomials. The practical problem is, however, that the matrix $`M^{(\alpha )}`$ is not well conditioned because it has eigenvalues with very different magnitudes. In order to illustrate this let us consider the special case $`(\alpha =1,\lambda =1,ϵ=0)`$ with $`n=10`$. In this case the eigenvalues are:

$`0.4435021205e14,0.1045947635e11,0.1143819915e9,`$   
$`0.7698917100e8,0.3571195735e6,0.1211873623e4,`$   
$`0.3120413130e3,0.6249495675e2,0.9849331094e1,`$   
$`1.075807246.`$ (8)   
A numerical investigation shows that, in general, the ratio of maximal to minimal eigenvalues is of the order of $`𝒪(10^{1.5n})`$. It is obvious from the structure of $`M^{(\alpha )}`$ in (5) that a rescaling of the interval $`[ϵ,\lambda ]`$ does not help. The large differences in magnitude of the eigenvalues implies through (6) large differences of magnitude in the coefficients $`c_{n\nu }`$ and therefore the numerical evaluation of the optimal polynomial $`P_n(x)`$ for large $`n`$ is non-trivial.

Let us now return to the general case with arbitrary function $`f(x)`$ and weight $`w(x)`$. It is very useful to introduce orthogonal polynomials $`\mathrm{\Phi }_\mu (x)(\mu =0,1,2,\mathrm{})`$ satisfying

$$_ϵ^\lambda 𝑑xw(x)^2\mathrm{\Phi }_\mu (x)\mathrm{\Phi }_\nu (x)=\delta _{\mu \nu }q_\nu .$$
(9)   
and expand the polynomial $`P_n(x)`$ in terms of them:

$$P_n(x)=\underset{\nu =0}{\overset{n}{}}d_{n\nu }\mathrm{\Phi }_\nu (x).$$
(10)   
Besides the normalization factor $`q_\nu `$ let us also introduce, for later purposes, the integrals $`p_\nu `$ and $`s_\nu `$ by

$`q_\nu `$ $``$ $`{\displaystyle _ϵ^\lambda }𝑑xw(x)^2\mathrm{\Phi }_\nu (x)^2,`$   
$`p_\nu `$ $``$ $`{\displaystyle _ϵ^\lambda }𝑑xw(x)^2\mathrm{\Phi }_\nu (x)^2x,`$   
$`s_\nu `$ $``$ $`{\displaystyle _ϵ^\lambda }𝑑xw(x)^2x^\nu .`$ (11)   
It can be easily shown that the expansion coefficients $`d_{n\nu }`$ minimizing $`\delta _n`$ are independent of $`n`$ and are given by

$$d_{n\nu }d_\nu =\frac{b_\nu }{q_\nu },$$
(12)   
where

$$b_\nu _ϵ^\lambda 𝑑xw(x)^2f(x)\mathrm{\Phi }_\nu (x).$$
(13)   
The minimal value of $`\delta _n^2`$ is

$$\delta _n^2=1N_{ϵ,\lambda }^1\underset{\nu =0}{\overset{n}{}}d_\nu b_\nu .$$
(14)   
Rescaling the variable $`x`$ by $`x^{}=\rho x`$ allows for considering only standard intervals, say $`[ϵ/\lambda ,1]`$. The scaling properties of the optimized polynomials can be easily obtained from the definitions. Let us now again consider the simple function $`f(x)=x^\alpha `$ and relative deviation with $`w(x)=x^\alpha `$ when the rescaling relations are:

$`\delta _n^2(\alpha ;ϵ\rho ,\lambda \rho )`$ $`=`$ $`\delta _n^2(\alpha ;ϵ,\lambda ),`$   
$`P_n(\alpha ;ϵ\rho ,\lambda \rho ;x)`$ $`=`$ $`\rho ^\alpha P_n(\alpha ;ϵ,\lambda ;x/\rho ),`$   
$`c_{n\nu }(\alpha ;ϵ\rho ,\lambda \rho )`$ $`=`$ $`\rho ^{\nu n\alpha }c_{n\nu }(\alpha ;ϵ,\lambda ).`$ (15)   
In applications to multi-bosonic algorithms for fermions the decomposition of the optimized polynomials as a product of root-factors is needed. This can be written as

$$P_n(\alpha ;ϵ,\lambda ;x)=c_{n0}(\alpha ;ϵ,\lambda )\underset{j=1}{\overset{n}{}}[xr_{nj}(\alpha ;ϵ,\lambda )].$$
(16)   
The rescaling properties here are:

$`c_{n0}(\alpha ;ϵ\rho ,\lambda \rho )`$ $`=`$ $`\rho ^{n\alpha }c_{n0}(\alpha ;ϵ,\lambda ),`$   
$`r_{nj}(\alpha ;ϵ\rho ,\lambda \rho )`$ $`=`$ $`\rho r_{nj}(\alpha ;ϵ,\lambda ).`$ (17)   
The root-factorized form (16) can also be used for the numerical evaluation of the polynomials with matrix arguments if a suitable optimization of the ordering of roots is performed .

The above orthogonal polynomials satisfy three-term recurrence relations which are very useful for numerical evaluation. In fact, at large $`n`$ the recursive evaluation of the polynomials is numerically more stable than the evaluation with root factors. For general $`f(x)`$ and $`w(x)`$, the first two ortogonal polynomials with $`\mu =0,1`$ are given by

$$\mathrm{\Phi }_0(x)=1,\mathrm{\Phi }_1(x)=x\frac{s_1}{s_0}.$$
(18)   
The higher order polynomials $`\mathrm{\Phi }_\mu (x)`$ for $`\mu =2,3,\mathrm{}`$ can be obtained from the recurrence relation

$$\mathrm{\Phi }_{\mu +1}(x)=(x+\beta _\mu )\mathrm{\Phi }_\mu (x)+\gamma _{\mu 1}\mathrm{\Phi }_{\mu 1}(x),(\mu =1,2,\mathrm{}),$$
(19)   
where the recurrence coefficients are given by

$$\beta _\mu =\frac{p_\mu }{q_\mu },\gamma _{\mu 1}=\frac{q_\mu }{q_{\mu 1}}.$$
(20)   
Defining the polynomial coefficients $`f_{\mu \nu }(0\nu \mu )`$ by

$$\mathrm{\Phi }_\mu (x)=\underset{\nu =0}{\overset{\mu }{}}f_{\mu \nu }x^{\mu \nu }$$
(21)   
the above recurrence relations imply the normalization convention

$$f_{\mu 0}=1,(\mu =0,1,2,\mathrm{}).$$
(22)   
The rescaling relations for the orthogonal polynomials easily follow from the definitions. For the simple function $`f(x)=x^\alpha `$ and relative deviation with $`w(x)=x^\alpha `$ we have

$$\mathrm{\Phi }_\mu (\alpha ;\rho ϵ,\rho \lambda ;x)=\rho ^\mu \mathrm{\Phi }_\mu (\alpha ;ϵ,\lambda ;x/\rho ).$$
(23)   
For the quantities introduced in (2) this implies

$`q_\nu (\alpha ;\rho ϵ,\rho \lambda )`$ $`=`$ $`\rho ^{2\alpha +1+2\nu }q_\nu (\alpha ;ϵ,\lambda ),`$   
$`p_\nu (\alpha ;\rho ϵ,\rho \lambda )`$ $`=`$ $`\rho ^{2\alpha +2+2\nu }p_\nu (\alpha ;ϵ,\lambda ),`$   
$`s_\nu (\alpha ;\rho ϵ,\rho \lambda )`$ $`=`$ $`\rho ^{2\alpha +1+\nu }s_\nu (\alpha ;ϵ,\lambda ).`$ (24)   
For the expansion coefficients defined in (12)-(13) one obtains

$`b_\nu (\alpha ;\rho ϵ,\rho \lambda )`$ $`=`$ $`\rho ^{\alpha +1+\nu }b_\nu (\alpha ;ϵ,\lambda ),`$   
$`d_\nu (\alpha ;\rho ϵ,\rho \lambda )`$ $`=`$ $`\rho ^{\alpha \nu }d_\nu (\alpha ;ϵ,\lambda ),`$ (25)   
and the recurrence coefficients in (19)-(20) satisfy

$`\beta _\mu (\alpha ;\rho ϵ,\rho \lambda )`$ $`=`$ $`\rho \beta _\mu (\alpha ;ϵ,\lambda ),`$   
$`\gamma _{\mu 1}(\alpha ;\rho ϵ,\rho \lambda )`$ $`=`$ $`\rho ^2\gamma _{\mu 1}(\alpha ;ϵ,\lambda ).`$ (26)   
For general intervals $`[ϵ,\lambda ]`$ and/or functions $`f(x)=x^\alpha \overline{P}(x)`$ the orthogonal polynomials and expansion coefficients have to be determined numerically. In some special cases, however, the polynomials can be related to some well know ones. An example is the weight factor

$$w^{(\rho ,\sigma )}(x)^2=(xϵ)^\rho (\lambda x)^\sigma .$$
(27)   
Taking, for instance, $`\rho =2\alpha ,\sigma =0`$ this weight is similar to the one for relative deviation from the function $`f(x)=x^\alpha `$, which would be just $`x^{2\alpha }`$. In fact, for $`ϵ=0`$ these are exactly the same and for small $`ϵ`$ the difference is negligible. The corresponding orthogonal polynomials are simply related to the Jacobi polynomials , namely

$$\mathrm{\Phi }_\nu ^{(\rho ,\sigma )}(x)=(\lambda ϵ)^\nu \nu !\frac{\mathrm{\Gamma }(\rho +\sigma +\nu +1)}{\mathrm{\Gamma }(\rho +\sigma +2\nu +1)}P_\nu ^{(\sigma ,\rho )}\left(\frac{2x\lambda ϵ}{\lambda ϵ}\right).$$
(28)   
Comparing different approximations with different $`(\rho ,\sigma )`$ the best choice is usually $`\rho =2\alpha ,\sigma =0`$ which corresponds to optimizing the relative deviation (see the appendix of ).

For large condition numbers $`\lambda /ϵ`$ least-squares optimization is much better than the Chebyshev approximation used for the approximation of $`x^1`$ in . The Chebyshev approximation is minimizing the maximum of the relative deviation

$$R(x)xP(x)1.$$
(29)   
For the deviation norm

$$\delta _{max}\underset{x[ϵ,\lambda ]}{\mathrm{max}}|R(x)|$$
(30)   
the least-squares approximation is slightly worse than the Chebyshev approximation. An example is shown by fig. 1. In the left lower corner the Chebyshev approximation has $`Rc(0.0002)=0.968`$ compared to $`Ro(0.0002)=0.991`$ for the least-squares optimization. For smaller condition numbers the Chebyshev approximation is not as bad as is shown by fig. 1. Nevertheless, in QCD simulations in sufficiently large volumes the condition number is of the order of the light quark mass squared in lattice units which can be as large as $`𝒪(10^610^7)`$.

Figure 1 also shows that the least-squares optimization is quite good in the minimax norm in (30), too. It can be proven that

$$|Ro(ϵ)|=\underset{x[ϵ,\lambda ]}{\mathrm{max}}|Ro(x)|$$
(31)   
hence the minimax norm can also be easily obtained from

$$\delta _{max}^{(o)}=\underset{x[ϵ,\lambda ]}{\mathrm{max}}|Ro(x)|=|Ro(ϵ)|.$$
(32)   
Therefore the least squares-optimization is also well suited for controlling the minimax norm, if for some reason it is required.

In QCD simulations the inverse power to be approximated ($`\alpha `$) is related to the number of Dirac fermion flavours: $`\alpha =N_f/2`$. If only $`u`$\- and $`d`$-quarks are considered we have $`N_f=2`$ and the function to be approximated is $`x^1`$. The dependence of the (squared) least-squares norm in (1) on the polynomial order $`n`$ is shown by fig. 2 for different values of the condition number $`\lambda /ϵ`$. The dependence on $`\alpha =N_f/2`$ is illustrated by fig. 3.

Another possible application of least-squares optimized polynomials is the numerical evaluation of the zero mass lattice action proposed by Neuberger . If one takes, for instance, the weight factor in (27) corresponding to the relative deviation, then the function $`x^{1/2}`$ has to be expanded in the Jacobi polynomials $`P^{(1,0)}`$.

## 3 Recurrence scheme

The expansion in orthogonal polynomials is very useful because it allows for a numerically stable evaluation of the least-squares optimized polynomials by the recurrence relation (19). The orthogonal polynomials themselves can also be determined recursively.

A recurrence scheme for obtaining the recurrence coefficients $`\beta _\mu ,\gamma _{\mu 1}`$ and expansion coefficients $`d_\nu =b_\nu /q_\nu `$ has been given in . In order to obtain $`q_\nu ,p_\nu `$ contained in (20) one can use the relations

$$q_\mu =\underset{\nu =0}{\overset{\mu }{}}f_{\mu \nu }s_{2\mu \nu },p_\mu =\underset{\nu =0}{\overset{\mu }{}}f_{\mu \nu }\left(s_{2\mu +1\nu }+f_{\mu 1}s_{2\mu \nu }\right).$$
(33)   
The coefficients themselves can be calculated from $`f_{11}=s_1/s_0`$ and (19) which gives

$`f_{\mu +1,1}`$ $`=`$ $`f_{\mu ,1}+\beta _\mu ,`$   
$`f_{\mu +1,2}`$ $`=`$ $`f_{\mu ,2}+\beta _\mu f_{\mu ,1}+\gamma _{\mu 1},`$   
$`f_{\mu +1,3}`$ $`=`$ $`f_{\mu ,3}+\beta _\mu f_{\mu ,2}+\gamma _{\mu 1}f_{\mu 1,1},`$   
$`\mathrm{}`$   
$`f_{\mu +1,\mu }`$ $`=`$ $`f_{\mu ,\mu }+\beta _\mu f_{\mu ,\mu 1}+\gamma _{\mu 1}f_{\mu 1,\mu 2},`$   
$`f_{\mu +1,\mu +1}`$ $`=`$ $`\beta _\mu f_{\mu ,\mu }+\gamma _{\mu 1}f_{\mu 1,\mu 1}.`$ (34)   
The orthogonal polynomial and recurrence coefficients are recursively determined by (20) and (33)-(3). The expansion coefficients for the optimized polynomial $`P_n(x)`$ can be obtained from

$$b_\mu =\underset{\nu =0}{\overset{\mu }{}}f_{\mu \nu }_ϵ^\lambda 𝑑xw(x)^2f(x)x^{\mu \nu }.$$
(35)   
The ingredients needed for this recursion are the basic integrals $`s_\nu `$ defined in (2) and

$$t_\nu _ϵ^\lambda 𝑑xw(x)^2f(x)x^\nu .$$
(36)   
The recurrence scheme based on the coefficients of the orthogonal polynomials $`f_{\mu \nu }`$ in (3) is not optimal for large orders $`n`$, neither for arithmetics nor for storage requirements. A better scheme can be built up on the basis of the integrals

$`r_{\mu \nu }{\displaystyle _ϵ^\lambda }𝑑xw(x)^2\mathrm{\Phi }_\mu (x)x^\nu ,`$   
$`\mu =0,1,\mathrm{},n;\nu =\mu ,\mu +1,\mathrm{},2n\mu .`$ (37)   
The recurrence coefficients $`\beta _\mu ,\gamma _{\mu 1}`$ can be expressed from

$$q_\mu =r_{\mu \mu },p_\mu =r_{\mu ,\mu +1}+f_{\mu 1}r_{\mu \mu }$$
(38)   
and eq. (20) as

$$\beta _\mu =f_{\mu 1}\frac{r_{\mu ,\mu +1}}{r_{\mu \mu }},\gamma _{\mu 1}=\frac{r_{\mu \mu }}{r_{\mu 1,\mu 1}}.$$
(39)   
It follows from the definition that

$`r_{0\nu }`$ $`=`$ $`{\displaystyle _ϵ^\lambda }𝑑xw(x)^2x^\nu =s_\nu ,`$   
$`r_{1\nu }`$ $`=`$ $`{\displaystyle _ϵ^\lambda }𝑑xw(x)^2(x^{\nu +1}+f_{11}x^\nu )=s_{\nu +1}+f_{11}s_\nu .`$ (40)   
The recurrence relation (19) for the orthogonal polynomials implies

$$r_{\mu +1,\nu }=r_{\mu ,\nu +1}+\beta _\mu r_{\mu \nu }+\gamma _{\mu 1}r_{\mu 1,\nu }.$$
(41)   
This has to be supplemented by

$$f_{11}=\frac{s_1}{s_0}$$
(42)   
and by the first equation from (3):

$$f_{\mu +1,1}=f_{\mu ,1}+\beta _\mu .$$
(43)   
Eqs. (39)-(43) define a complete recurrence scheme for determining the orthogonal polynomials $`\mathrm{\Phi }_\mu (x)`$. The moments $`s_\nu `$ of the integration measure defined in (2) serve as the basic input in this scheme.

The integrals $`b_\nu `$ in (13), which are necessary for the expansion coefficients $`d_\nu `$ in (12), can also be calculated in a similar scheme built up on the integrals

$`b_{\mu \nu }{\displaystyle _ϵ^\lambda }𝑑xw(x)^2f(x)\mathrm{\Phi }_\mu (x)x^\nu ,`$   
$`\mu =0,1,\mathrm{},n;\nu =0,1,\mathrm{},n\mu .`$ (44)   
The relations corresponding to (3)-(41) are now

$`b_{0\nu }={\displaystyle _ϵ^\lambda }𝑑xw(x)^2f(x)x^\nu =t_\nu ,`$ (45)   
$`b_{1\nu }={\displaystyle _ϵ^\lambda }𝑑xw(x)^2f(x)(x^{\nu +1}+f_{11}x^\nu )=t_{\nu +1}+f_{11}t_\nu ,`$   
$`b_{\mu +1,\nu }=b_{\mu ,\nu +1}+\beta _\mu b_{\mu \nu }+\gamma _{\mu 1}b_{\mu 1,\nu }.`$   
The only difference compared to (3)-(41) is that the moments of $`w(x)^2`$ are now replaced by the ones of $`w(x)^2f(x)`$.

It is interesting to collect the quantities which have to be stored in order that the recurrence can be resumed. This is useful if after stopping the iterations, for some reason, the recurrence will be restarted. Let us assume that the quantities $`q_\nu (\nu =0,\mathrm{},n)`$, $`b_\nu (\nu =0,\mathrm{},n)`$, $`\beta _\nu (\nu =1,\mathrm{},n1)`$ and $`\gamma _\nu (\nu =0,\mathrm{},n2)`$ are already known and one wants to resume the recurrence in order to calculate these quantities for higher indices. For this it is enough to know the values of

$`f_{n1,1},r_{n1,n1},`$   
$`R_{0\mathrm{}n}^{(0)}(r_{0,2n+1},r_{1,2n},\mathrm{},r_{n,n+1}),`$   
$`R_{0\mathrm{}n}^{(1)}(r_{0,2n},r_{1,2n1},\mathrm{},r_{n,n}),`$   
$`B_{0\mathrm{}n}^{(1)}(b_{0,n},b_{1,n1},\mathrm{},b_{n,0}),`$   
$`B_{0\mathrm{}n1}^{(2)}(b_{0,n1},b_{1,n2},\mathrm{},b_{n1,0}).`$ (46)   
This shows that for maintaining a resumable recurrence it is enough to store a set of quantities linearly increasing in $`n`$.

An interesting question is the increase of computational load as a function of the highest required order $`n`$. At the first sight this seems to go just like $`n^2`$, which is surprising because, as eq. (6) shows, finding the minimum requires the inversion of an $`nn`$ matrix. However, numerical experience shows that the number of required digits for obtaining a precise result does also increase linearly with $`n`$. This is due to the linearly increasing logarithmic range of eigenvalues, as illustrated by (2). Using, for instance, Maple V for the arbitrary precision arithmetic, the computation slows down by another factor going roughly as (but somewhat slower than) $`n^2`$. Therefore, the total slowing down in $`n`$ is proportional to $`n^4`$. For the same reason the storage requirements increase by $`n^2`$.

## 4 A convenient choice for TSMB

In the TSMB algorithm for Monte Carlo simulations of fermionic theories, besides the simple function $`x^\alpha `$, also the function $`x^\alpha /\overline{P}(x)`$ has to be approximated. Here $`\overline{P}(x)`$ is typically a lower order approximation to $`x^\alpha `$. In this case, if one chooses to optimize the relative deviation, the basic integrals defined in (2) and (36) are, respectively,

$`s_\nu `$ $`=`$ $`{\displaystyle _ϵ^\lambda }𝑑x\overline{P}(x)^2x^{2\alpha +\nu },`$   
$`t_\nu `$ $`=`$ $`{\displaystyle _ϵ^\lambda }𝑑x\overline{P}(x)x^{\alpha +\nu }.`$ (47)   
It is obvious that, if the recurrence coefficients for the expansion of the polynomial $`\overline{P}(x)`$ in orthogonal polynomials are known, the recursion scheme can also be used for the evaluation of $`s_\nu `$ and $`t_\nu `$.

Another observation is that the integrals in (4) can be simplifyed if, instead of choosing the weight factor $`w(x)^2=\overline{P}(x)^2x^{2\alpha }`$, one takes

$$w(x)^2=\overline{P}(x)x^\alpha ,$$
(48)   
which leads to

$`s_\nu `$ $`=`$ $`{\displaystyle _ϵ^\lambda }𝑑x\overline{P}(x)x^{\alpha +\nu },`$   
$`t_\nu `$ $`=`$ $`{\displaystyle _ϵ^\lambda }𝑑xx^\nu .`$ (49)   
Since $`\overline{P}(x)`$ is an approximation to $`x^\alpha `$, the function $`f(x)x^\alpha /\overline{P}(x)`$ is close to one and the difference between $`f(x)^2`$ and $`f(x)^1`$ is small. Therefore the least-squares optimized approximations with the weights $`w(x)^2=f(x)^2`$ and $`w(x)^2=f(x)^1`$ are also similar. It turns out that the second choice is, in fact, a little bit better because the largest deviation from $`x^\alpha `$ typically occurs at the lower end of the interval $`x=ϵ`$ where $`\overline{P}(x)`$ is smaller than $`x^\alpha `$. As a consequence, $`\overline{P}(x)x^\alpha <1`$ and $`\overline{P}(x)x^\alpha >\overline{P}(x)^2x^{2\alpha }`$. This means that choosing the weight factor in (48) is emphasising more the lower end of the interval where $`\overline{P}(x)`$ as an approximation of $`x^\alpha `$ is worst.

In summary: least-squares optimization is a flexible and powerful tool which can serve as a basis for applying the two-step multi-bosonic algorithm for Monte Carlo simulations of QCD and other similar theories. With the help of the recurrence scheme described in the previous section one can determine the necessary polynomial approximations to high enough orders.