=Paper=
{{Paper
|id=Vol-3754/paper10
|storemode=property
|title=Solving estimation problems using minimax polynomials and Gröbner bases
|pdfUrl=https://ceur-ws.org/Vol-3754/paper10.pdf
|volume=Vol-3754
|authors=Kenta Kuramochi,Akira Terui,Masahiko Mikawa
|dblpUrl=https://dblp.org/rec/conf/sycss/KuramochiTM24
}}
==Solving estimation problems using minimax polynomials and Gröbner bases==
https://ceur-ws.org/Vol-3754/paper10.pdf
Solving Estimation Problems Using Minimax Polynomials
and Gröbner Bases⋆
Kenta Kuramochi1 , Akira Terui2,∗ and Masahiko Mikawa3
1
Master’s Program in Mathematics, Graduate School of Science and Technology, University of Tsukuba, Tsukuba 305-8571, Japan
2
Institute of Pure and Applied Sciences, University of Tsukuba, Tsukuba 305-8571, Japan
3
Institute of Library, Information and Media Science, University of Tsukuba, Tsukuba 305-8550, Japan
Abstract
We propose a method for solving the speech direction estimation problem by computer algebra. The method is
based on the function approximation using the minimax polynomial. The minimax polynomial is obtained by an
iterative method called the Remez exchange algorithm, in which Gröbner bases computation is employed. We
present an effective way to compute the minimax polynomial using Gröbner bases.
Keywords
Estimation problem, Function approximation, Minimax polynomial, Remez exchange algorithm, Gröbner basis,
Speech direction estimation problem
1. Intorduction
In this paper, we discuss solving estimation problems with a function approximation method using the
minimax polynomial and computer algebra.
Function approximation is the technique to approximate functions. It is used to make a sequence of
polynomials for proving the density of function space [1] or to regard a function as a polynomial for
evaluation [2]. Various methods for function approximation have been proposed, such as the Maclaulin
expansion or the least squares method [2]. Here, as one of them, we present the minimax approximation
and Remez exchange algorithm. The minimax approximation is the approximation using the polynomial
which meets the property that the maximum value of the difference between the given function and
the derived polynomial is the smallest of all polynomials in a given domain. The polynomial satisfying
such a property is called the minimax polynomial. Since minimax polynomials are polynomials, one
can use algebraic computation. On the other hand, the interpolation method [2], which is frequently
used in computer algebra, estimates the polynomial based on discrete points and values on the original
function. Compared with the interpolation method, the minimax approximation is better for errors, that
is, the maximum value of errors of the minimax polynomial is less than that of the polynomial obtained
by the interpolation method in many cases. In computing the minimax polynomial, an iteration method
called the Remez exchange algorithm is used.
Estimation problems are the problems of estimating unknown information using already known
information. Solving estimation problems is important in developing devices since objects that are
measurable are limited. To solve estimation problems, one uses numerical methods such as the gradient
method [3] or the genetic Algorithms [4]. However, the gradient method may return a local solution
depending on initial values since it uses local convergence properties, and the genetic algorithm has some
disadvantages in that it sometimes solves the estimation problem not properly due to the phenomena
called initial convergence and hitchhiking. On the other hand, the estimation method using minimax
SCSS 2024: 10th International Symposium on Symbolic Computation in Software Science, August 28–30, 2024, Tokyo, Japan
⋆
This work was partially supported by JKA and its promotion funds from KEIRIN RACE.
∗
Corresponding author.
Envelope-Open s2320136@u.tsukuba.ac.jp (K. Kuramochi); terui@math.tsukuba.ac.jp (A. Terui); mikawa@slis.tsukuba.ac.jp (M. Mikawa)
GLOBE https://researchmap.jp/aterui (A. Terui); https://mikawalab.org/ (M. Mikawa)
Orcid 0000-0003-0846-3643 (A. Terui); 0000-0002-2193-3198 (M. Mikawa)
© 2024 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
CEUR
ceur-ws.org
Workshop ISSN 1613-0073
Proceedings
57
approximation together with Gröbner bases [5] computation may avoid these phenomena, for this
method evaluates values globally.
In this paper, we propose a method for solving estimation problems using the minimax polynomial
and Gröbner bases, as follows. First, for a given mathematical model (function) of the estimation
problem, we calculate the minimax polynomial which approximates the given function. Then, we
make a system of polynomial equations and solve it with Gröbner bases computation for obtaining
the solution of the estimation problem. Furthermore, we apply the proposed method to the speech
direction estimation problem, which is an estimation of the direction of a speaker using a microphone
array. For the speech direction estimation problem, a method using the Genetic Algorithm has been
proposed [6]. We show that we can effectively use the minimax approximation and Gröbner bases for
finding global solutions to the speech direction estimation problem.
The paper is organized as follows. In Section 2, we review the definition of the minimax polynomial
and the Remez exchange algorithm. In Section 3, after defining the estimation problem more minutely,
we propose a solution for this using the minimax approximation and Gröbner bases. In Section 4, we
introduce the speech direction estimation problem, which is the task we are working on, and explain
the reason why the method in this paper can be used. In Section 5, we conclude and pick up some
challenges we are facing and tasks to improve our tasks.
2. Preliminaries
Let 𝐾 be a field. The notation 𝐾 [𝑥] or 𝐾 [𝑥1 , … , 𝑥𝑛 ] stands for the ring of polynomials over 𝐾 in 𝑥1 , … , 𝑥𝑛 .
For a function 𝑓, ‖𝑓 ‖∞ denotes the infinity norm of 𝑓.
2.1. Minimax Approximation
In the following, let [𝑎, 𝑏] be a closed interval and 𝑓 be a continuous function on [𝑎, 𝑏].
Definition 2.1 (Minimax polynomial). For the function 𝑓 in above, a polynomial 𝑃 ∈ 𝐾 [𝑥] of degree 𝑘
which minimizes
‖𝑓 − 𝑃‖∞ = max |𝑓 (𝑥) − 𝑃(𝑥)| ,
𝑎≤𝑥≤𝑏
is called the 𝑘-th minimax polynomial of 𝑓.
Definition 2.1 says that, if 𝑃 ∈ 𝐾 [𝑥] is the 𝑘-th minimax polynomial of 𝑓, the inequality
max |𝑓 (𝑥) − 𝑃(𝑥)| ≤ max |𝑓 (𝑥) − 𝑄(𝑥)|,
𝑎≤𝑥≤𝑏 𝑎≤𝑥≤𝑏
follows for any 𝑄 ∈ 𝐾 [𝑥] with deg 𝑄 = 𝑘. In other words, the 𝑘-th minimax polynomial is the best
polynomial in terms of error. Thus, one should use the minimax polynomial if one wants to approximate
the given function by polynomials with minimizing errors.
The following theorem [2] tells us that that that value has a minimum value.
Theorem 2.1. If the function 𝑓 is continuous, the 𝑓 has a unique 𝑘-th minimax polynomial for any 𝑘 ∈ ℤ≥0 .
Furthermore, we can construct the 𝑘-th degree minimax polynomial by Algorithm 1.
For details of Algorithm 1 such as uniqueness and termination, see [7]. Note that Algorithm 1
needs to compute the solution of a system of linear equations and the extremum points of continuous
functions. (We have implemented Algorithm 1 using a computer algebra system Risa/Asir [8] with the
library os_muldif [9] 1 .)
The method to approximate functions using the Remez exchange algorithm is called minimax
approximation. Note that the maximum error of 𝑓 − 𝑃 becomes smaller as the degree of minimax
1
Solving a system of linear equations is performed by Risa/Asir itself and the function”os_md.fmmx()” in the os_muldif library
is used for computing extremum points with the extrema.
58
Algorithm 1 Remez exchange algorithm [2]
Require: a continuous function 𝑓, a closed interval [𝑎, 𝑏], a degree of a polynomial𝑘
Ensure: a minimax polynomial 𝑃 ∈ 𝐾 [𝑥] with degree 𝑘
1: 𝐸max ≔ ∞, 𝐸min ≔ 0
2: 𝐼 ≔ [𝑥0 , 𝑥1 , … , 𝑥𝑘 , 𝑥𝑘+1 ] : 𝑘 + 2 initial interpolate points
3: while 𝐸max >> 𝐸min do
𝑘
4: [𝑎0 , 𝑎1 , … , 𝑎𝑘 , 𝐸] : the solution of ∑𝑗=0 𝑎𝑗 (𝑥𝑖 )𝑗 + (−1)𝑖 𝐸 = 𝑓 (𝑥𝑖 ) (𝑖 = 0, 1, … , 𝑘, 𝑘 + 1)
𝑘
5: 𝑃 ≔ ∑𝑗=0 𝑎𝑗 𝑥 𝑗
6: 𝐼 ≔ [𝑥0 , 𝑥1 , … , 𝑥𝑘 , 𝑥𝑘+1 ] : the 𝑘 + 2 extremum points of 𝐸(𝑥) ≔ 𝑓 (𝑥) − 𝑃(𝑥)
7: 𝐸max ≔ max{|𝐸(𝑥0 )|, |𝐸(𝑥1 )|, … , |𝐸(𝑥𝑘+1 )|}, 𝐸min ≔ min{|𝐸(𝑥0 )|, |𝐸(𝑥1 )|, … , |𝐸(𝑥𝑘+1 )|}
8: end while
9: return 𝑃
polynomial 𝑘 becomes larger and the interval [𝑎, 𝑏] more restricted. The method in this paper can
be used if functions are continuous, their range is bounded and the variable to solve is in a bounded
interval. If the variable to solve is in a bounded interval that is not closed, we can construct the minimax
polynomial over a closed one containing it.
Note that the minimax polynomial 𝑃 of 𝑓 has the following property: the error function 𝑓 − 𝑃 has
deg 𝑃 + 2 maximum value with alternate change of singers in [𝑎, 𝑏]. The value of |𝑓 (𝑥) − 𝑃(𝑥)| rapidly
becomes larger as the value 𝑥 separates from [𝑎, 𝑏].
3. The estimation problem
In this section, we introduce the method to solve estimation problems using minimax approximation
and computer algebra.
Let 𝑢 = (𝑢1 , … , 𝑢𝑛 ) be measurable parameters and 𝑣 = (𝑣1 , … , 𝑣𝑡 ) be immeasurable parameters, and 𝑥
be the parameter to be estimated. Assume that a mathematical model describing phenomena is given as
𝑓 (𝑢, 𝑥) + ∑ 𝑔𝑖 (𝑢)ℎ𝑗 (𝑣𝑘 ) = 0, (1)
𝑖,𝑗,𝑘
where 𝑓 is a continuous function and 𝑥 is bounded. Since 𝑢 are measurable, one can consider 𝑓 (𝑢, 𝑥) as
a function only in 𝑥, 𝑔𝑖 (𝑢) as a constant by substituting 𝑢 and ℎ𝑖 (𝑣𝑗 ) as a new immeasurable parameter
by replacing it with a new variable properly if necessary. Thus, 𝑓 can be approximated by the minimax
polynomial and the equation is transformed into the form of
𝑃(𝑥) + ∑ 𝑄𝑖 (𝑣) = 0, (2)
𝑖
where 𝑃 and 𝑄𝑖 represents the polynomials in 𝑥 and 𝑣, respectively. As eq. (2) is a multivariate polynomial
equation, a system of polynomial equations can be generated by substituting measured values into
measurable parameters. Thus, Gröbner bases and the Elimination Theorem [5] can be used to solve the
system.
Note that the minimax approximation is the best way to approximate a function by polynomials in
terms of errors. However, under some conditions, the minimax approximation cannot be performed. To
be able to apply the approximation, the mathematical model with the estimated parameter must be
continuous, have no other immeasurable parameters, and the estimated parameter must be in a bounded
interval. The more restricted the interval is, the better the accuracies of the minimax approximations are.
Thus this method is suitable for estimation problems in that the range of the solution of the estimated
parameters is bounded.
59
4. An example by speech direction estimation problem
Speech direction estimation problem [6] is the problem of estimating the orientation of the speaker’s
face. This problem derives from robotics. Nowadays robots have been playing important roles in various
fields and some research on robots cooperating with each other is underway [10, 11]. For such tasks,
how to select one robot from more than one is one of the themes when humans select and command
one to do some tasks. Thus, we are trying to select one robot by calling like humans saying “excuse
me” or “hey.” We assume a situation in a room like we are calling one robot waiter from many in a
restaurant. Thus naming each robot is difficult.
We suppose that we know the coordinates of a user and robots, for the microphone arrays we are
assuming have a function to compute the locations of themselves and the directions of arrival (DoA)
and we can compute the user’s location using DoA. For details, see [6]. To estimate the direction of the
speaker’s face orientation, we use a mathematical model called the voice spread model. The voice spread
model is a formula to compute logical sound pressure levels recorded by a particular microphone array.
We consider the situation in which there is a speaker and microphone array mic 𝑖. To construct a voice
spread model, we need to consider two attenuation effects, distance attenuation and angle attenuation.
Distance attenuation is the effect that the longer the distance between a sound source and a sound
receiving point is, the smaller the sound pressure level of the receiving point is. Supposing a sound
source to be a point, the sound pressure level is in inverse proportion to the square of the distance.
Angle attenuation is the effect that the sound pressure level is the strongest in the front direction of
a mouth and it gets weaker when separate from the direction. There are various kinds of research on
angle attenuations [12, 13, 14] and we adopt Monson’s [15].
Let Σ𝑊 be the world coordinate system, Σ𝑆 be the local coordinate system whose center is the user and
the 𝑥-axis is the front direction of the user’s face. Furthermore, let 𝜃 be the angle formed by the 𝑥-axis
of Σ𝑊 and Σ𝑆 . Note that 𝜃 is the very estimated parameter. Then, the coordinate of mic 𝑖 is expressed as
𝑆 𝑝 = ( cos 𝜃 sin 𝜃
) (𝑊 𝑝𝑖 − 𝑊 𝑝𝑆 ).
𝑖 − sin 𝜃 cos 𝜃
The parameters 𝑊 𝑝𝑖 and 𝑊 𝑝𝑆 stand for the coordinate of mic 𝑖 and the user in Σ𝑊 , respectively. Then,
the theoretical sound pressure level 𝐿̂ 𝑖 recorded by mic 𝑖 is described as
1 + cos(𝜑𝑖 )
𝐿̂ 𝑖 = 𝐿 − 10 log10 ‖𝑆 𝑝𝑖 ‖2 − 𝑎 (1 − ).
2
Note that parameter 𝐿 is the sound pressure level of the point that is 1 [m] away from the speaker and in
the direction of the speaker’s mouth. A parameter 𝑎 is the size of angle attenuation and 𝜙𝑖 represents the
azimuth angle to mic 𝑖 in Σ𝑆 . Note that the parameters 𝐿 and 𝑎 are immeasurable and 𝜙𝑖 is measurable.
Let Σ𝑊 be the world coordinate system, Σ𝑆 be the local coordinate system whose center is the user
and 𝑥-axis is the front direction of the user’s face. Also let 𝜃 be the angle formed by 𝑥-axises of Σ𝑊 and
Σ𝑆 . Note that 𝜃 is the very estimated parameter. The coordinate of mic 𝑖’s microphone array is expressed
as
𝑆 𝑝 = ( cos 𝜃 sin 𝜃
) ( 𝑊 𝑝𝑖 − 𝑊 𝑝𝑆 )
𝑖 − sin 𝜃 cos 𝜃
The parameters 𝑊 𝑝𝑖 and 𝑊 𝑝𝑆 stands for the coordinate of mic 𝑖 and the user in Σ𝑊 respectively. Then,
the theoritical sound pressure level 𝐿̂ 𝑖 recorded by mic 𝑖 is described as
1 + cos(𝜑𝑖 )
𝐿̂ 𝑖 = 𝐿 − 10 log10 ‖𝑆 𝑝𝑖 ‖2 − 𝑎 (1 − )
2
Note that parameter 𝐿 is the sound pressure level of the point that is 1 [m] away from the speaker and
in the direction of the speaker’s mouse. A parameter 𝑎 is the size of angle attenuation and 𝜙𝑖 represents
the azimuth angle to mic 𝑖 in Σ𝑆 . Make sure that the parameters 𝐿 and 𝑎 are immeasurable and 𝜙𝑖 is
measurable.
60
Since 𝑊 𝑝𝑖 , 𝑊 𝑝𝑆 , and 𝜑𝑖 are measurable parameters, one can substitute measured values into these.
1 + cos(𝜑𝑖 )
Thus, the term (1 − ) can be alternated by a measured parameter 𝑚𝑖 and 𝐿̂ 𝑖 can be substituted
2
by a measured value. Thus, the equation above can be transformed into the form of
𝐿̂ 𝑖 − 𝐿 + 10 log10 ‖𝑆 𝑝𝑖 ‖2 + 𝑎𝑚𝑖 = 0.
Since the term log10 ‖𝑆 𝑝𝑖 ‖2 can be regarded as a function with the variable 𝜃 only and 𝜃 is in a bounded
interval (−𝜋, 𝜋], we can approximate that term by the minimax polynomial 𝑃𝑖 (𝜃). As a result, we have
𝐿̂ 𝑖 − 𝐿 + 10𝑃𝑖 (𝜃) + 𝑎𝑚𝑖 = 0.
The equation above is a multivariate polynomial in 𝐿, 𝜃, 𝑎 for any number 𝑖. To find 𝜃 (−𝜋 < 𝜃 ≤ 𝜋), we
need to measure parameters 𝑊 𝑝𝑖 , 𝑊 𝑝𝑆 , 𝜑𝑖 , 𝐿̂ 𝑖 for given 𝑖 ≥ 3, thus we compute the Gröbner bases and
make the system of equations a triangular form.
5. Concluding remarks
In this paper, we have proposed a method for estimation problems using the minimax approximation
and Gröbner bases computation. As an application, we have shown that the method can be used for the
speech direction estimation problem.
To check the accuracy of the proposed method, experiments should be carried out in simulations
and actual environments. For simulations of speech direction estimation problems, Python language
has a library called Pyroomacoustics [16]. It is a software library to develop and test algorithms for
voice processing, and can simulate an environment in a room and one can adjust the properties of the
environment such as the location and the directivity of a sound source and a microphone, temperature,
humidity, dimension, size and shape of a room and so on.
Furthermore, there exist the minimax rational functions and algorithms for computing ones. The
minimax rational function of 𝑓 is a rational function 𝑅𝑚,𝑘 which minimizes
𝑚
∑𝑗=0 𝑎𝑗 𝑥 𝑗
max |𝑓 (𝑥) − 𝑅𝑚,𝑘 (𝑥)| = max |𝑓 (𝑥) − |,
𝑎≤𝑥≤𝑏 𝑎≤𝑥≤𝑏 𝑘
∑𝑗=0 𝑏𝑗 𝑥 𝑗
for given 𝑚, 𝑘 ∈ ℤ≥1 . Compared with the minimax polynomials, the minimax rational functions tend
to have less errors. However, we failed to compute one properly due to errors in the calculation of
improper integrals. We need to compute Chebyshev expansion [2] of 𝑓 to construct a minimax rational
function and to compute improper integrals. In our experiment, we failed to compute improper integrals
properly due to divergence of integrands. Thus, the method to compute improper integrals with fewer
errors should be investigated.
To solve estimation problems independent of the values of measured parameters, we need to ap-
proximate multivariate functions with multivariate polynomials or rational functions. Cody [17] says
almost no algorithm for approximating those with minimax polynomials (or rational functions) exists.
However, Luke [18] has made approximations of a variety of functions in mathematical physics using
hypergeometric functions. Loeb [19] reports a linear algorithm to create rational approximations over
a discrete point set. Fox, Goldstein and Lastman [20] also have proposed an algorithm for rational
approximation on finite point sets.
The problem with which this method can be used is very limited since mathematical models that
satisfy equation (1) and estimated parameters to be bounded are very restrictive. However, direction
estimation problems are good because parameters of direction are usually bounded on (−𝜋, 𝜋] or [0, 2𝜋).
To see the efficacy or effectiveness of this method, we should seek estimation problems that can be
solved by this method.
61
References
[1] H. M. Stone, The generalized weierstrass approximation theorem, Mathematics Magazine 21
(1948) 167–184. doi:10.2307/3029750 .
[2] H. M. Antia, Numerical methods for scientists and engineers, second ed., Basel; Boston: Birkhäuser,
2002. doi:10.1007/978- 93- 86279- 52- 1 .
[3] E. Polak, Optimization: Algorithms and Consistent Approximations, Springer, 1997. doi:10.1007/
978- 1- 4612- 0663- 7 .
[4] M. Mitchell, An Introduction to Genetic Algorithms, The MIT Press, 1998. doi:10.7551/mitpress/
3927.001.0001 .
[5] D. Cox, J. Little, D. O’shea, M. Sweedler, Ideals, varieties, and algorithms: An Intorduction to
Computational Algebraic Geometry and Commutative Algebra, fourth ed., Springer, 2015. doi:10.
1007/978- 3- 319- 16721- 3 .
[6] T. Kato, Individual Selection Method for Multi-Robot based on Speech Direction Estimation (in
Japanese), Master’s thesis, Institute of Library, Information and Media Science, University of
Tsukuba, 2024. 82 pages.
[7] N. Daili, A. Guesmia, Remez algorithm applied to the best uniform polynomial approximations,
General Mathematics Notes 17 (2013) 16–31.
[8] M. Noro, T. Takeshima, Risa/Asir — A Computer Algebra System, in: ISSAC ’92: Papers from the
International Symposium on Symbolic and Algebraic Computation, Association for Computing
Machinery, New York, NY, USA, 1992, pp. 387–396. doi:10.1145/143242.143362 .
[9] T. Oshima, os_muldif: a library for computer algebra Risa/Asir, 2007–2024. URL: https://www.ms.
u-tokyo.ac.jp/~oshima/muldif/os_muldifeg.pdf, accessed 2024-06-05.
[10] P. R. Wurman, R. D’Andrea, M. Mountz, Coordinating hundreds of cooperative, autonomous
vehicles in warehouses, AI Magazine 29 (2008) 9–19. doi:10.1609/aimag.v29i1.2082 .
[11] J. Alonso-Mora, A. Breitenmoser, M. Rufli, R. Siegwart, P. Beardsley, Image and animation display
with multiple mobile robots, The International Journal of Robotics Research 132 (2012) 753–773.
doi:10.1177/0278364912442095 .
[12] W. T. Chu, A. Warnock, Detailed directivity of sound fields around human talkers, Technical
Report RR-104, National Research Council of Canada, 2002. doi:10.4224/20378930 .
[13] G. A. Studebaker, Directivity of the human vocal source in the horizontal plane, Ear and hearing 6
(1985) 315–319. doi:10.1097/00003446- 198511000- 00007 .
[14] D. Cebrera, P. J. Davis, A. Connolly, Long-term horizontal vocal directivity of opera singers:
Effects of singing projection and acoustic environment, Journal of Voice 25 (2011) e291–e303.
doi:10.1016/j.jvoice.2010.03.001 .
[15] B. B. Monson, E. J. Hunter, B. H. Story, Horizontal directivity of low-and high-frequency energy
in speech and singing, The Journal of the Acoustical Society of America 132 (2012) 433–411.
doi:10.1121/1.4725963 .
[16] R. Scheibler, E. Bezzam, I. Dokmanić, Pyroomacoustics: A python package for audio room
simulation and array processing algorithms, in: 2018 IEEE International Conference on Acoustics,
Speech and Signal Processing (ICASSP), IEEE, 2018, pp. 351–355. doi:10.1109/ICASSP.2018.
8461310 .
[17] W. J. Cody, A survey of practical rational and polynomial approximation of functions, SIAM
Review 12 (1970) 400–423. doi:10.1137/1012082 .
[18] Y. L. Luke, The special functions and their approximations, volume 1, Academic Press, 1969.
[19] H. L. Loeb, Algorithms for chebyshev approximation using the ratio of linear forms, Journal of
the Society for Industrial and Applied Mathematics 8 (1960) 458–465. doi:10.1137/0108031 .
[20] P. Fox, A. Goldstein, G. Lastman, Rational approximations on finite point sets, in: Approximation
of Functions (Proc. Sympos. General Motors Res. Lab., 1964), 1965, p. 57.
62