11institutetext: University of Coimbra, Centre for Informatics and Systems of
the University of Coimbra, Department of Informatics Engineering 11email:
<EMAIL_ADDRESS>
# Context Matters: Adaptive Mutation for Grammars

Pedro Carvalho  0000-0003-3845-4617 Jessica Mégane  0000-0001-6697-5423 Nuno
Lourenço  0000-0002-2154-0642 Penousal Machado  0000-0002-6308-6484

###### Abstract

This work proposes Adaptive Facilitated Mutation, a self-adaptive mutation
method for Structured Grammatical Evolution (SGE), biologically inspired by
the theory of facilitated variation. In SGE, the genotype of individuals
contains a list for each non-terminal of the grammar that defines the search
space. In our proposed mutation, each individual contains an array with a
different, self-adaptive mutation rate for each non-terminal. We also propose
Function Grouped Grammars, a grammar design procedure to enhance the benefits
of the propose mutation. Experiments were conducted on three symbolic
regression benchmarks using Probabilistic Structured Grammatical Evolution
(PSGE), a variant of SGE. Results show our approach is similar or better when
compared with the standard grammar and mutation.

###### Keywords:

Adaptive Mutation Grammar-design Grammar-based Genetic Programming.

##  1 Introduction

Grammar-based Genetic Programming (GP) algorithms have been an important tool
for the evolution of computer programs since their inception [1, 2, 3]. The
most popular approach is Grammatical Evolution (GE) which is notable for
decoupling the genotype and the phenotype, using a grammar to translate a data
structure into an executable program. The representation and variation
operators used by GE present some known issues, such as low locality and high
redundancy. The first means that small changes in the genotype can cause
significant changes in the phenotype, and the second means that most
modifications do not affect the phenotype. These characteristics result in a
bad trade between exploration and exploitation, which makes the algorithm
perform similarly to random search [4].

Structured Grammatical Evolution (SGE) [5] is a variant of GE that uses a
different representation for the individuals. The genotype comprises several
lists, one for each non-terminal in the grammar, and each list contains the
indexes of the rules to be expanded. SGE shows better performance when
compared to GE, and other grammar-based approaches [6], but also improved
locality and lower redundancy when compared to standard GE [7, 8], in part due
to its operators. This representation allows the recombination operator to be
grammar-aware, preserving the list of each non-terminal. On the other hand,
the grammar does not inform the mutation operator. Mutation in SGE changes the
production rule of the non-terminal selected to mutate. This operator affects
all genes with the same frequency, regardless of grammatical context. Using a
static and equal value for all non-terminals fails to consider that not all
mutations are equally destructive.

Specific genes can play an essential role in the solution and, when mutated,
may completely ruin the phenotype behavior. On the other hand, some genes may
have a tuning role; in this case, a mutation will only result in a minor
adjustment to the solution. Despite these differences, both types of genes are
equal in the eyes of mutation. Biological processes have evolved to prevent
this phenomenon. Gerhart et al. [9] propose that there are core components
vital to the individual which remain unchanged for long periods and regulatory
genes that combine existing core components and change frequently. The result
is a system that can quickly adapt to new environments through regulatory
changes while preserving the core components that ensure individuals are
functional. We can replicate this behavior in grammar-based algorithms using
different adaptive mutation probabilities for each non-terminal. This approach
enables the system to autonomously regulate mutation rates to match the impact
of changes to that non-terminal. Note that this solution is only as effective
as the correlation between non-terminals and mutation impact. It follows that
the effectiveness of this mutation is related to grammar-design [10, 11, 12,
13], as grammars with more rules enable finer tuning of mutation
probabilities. The grouping of productions within each non-terminal is also
relevant; separating low and high impact changes into separate rules should
improve performance.

In this work, we propose Adaptive Facilitated Mutation, a biologically
inspired grammar-aware self-adaptive mutation operator for SGE and its
variants. Furthermore, we propose ”Function Grouped Grammars”, a method for
grammar design that empirically outperforms grammars commonly used for
regression in GE. We compare our approach to standard grammar and mutation and
find that, when combined, Function Grouped Grammars and Adaptive Facilitated
Mutation are statistically superior or similar to the baseline in three
relevant GP benchmarks.

The remainder of this work is structured as follows: First, section 2 presents
the background necessary to understand the work presented. Section 3 presents
the proposed mutation and grammar-design method. Section 4 details the
experimentation setup used, and Section 5 the experimental results regarding
performance and analysis of probabilities. Section 6 gathers the main
conclusions and provides insights regarding future work.

##  2 Background

Evolutionary Algorithms (EAs) are optimization algorithms inspired by the
biological processes of natural evolution. A population of individuals
(candidate solutions) evolves over several generations, guided by a fitness
function. Similar to nature, these individuals are subject to selection,
reproduction, and genetic variation. These algorithms face known issues such
as parameter tuning, premature convergence, and lack of diversity. Researchers
propose novel representations [14, 5, 15], selection methods, genetic
operators, and parameter selection approaches to address these issues.

###  2.1 Grammar-based Genetic Programming

GP is an EA that evolves solutions as programs. Over the years, researchers
have proposed many variants of GP, and grammar-based approaches gained more
popularity as grammars are helpful to set restrictions to the search space
[3].

GE [14] is the most popular grammar-based GP methods. The individuals’
genotype is a string/vector of integers that is translated into a phenotype
(an executable function) through a grammar. The individuals are subject to
selection mechanisms, mutation, and crossover in each generation.

This approach is relevant but suffers from high redundancy [16], and poor
locality [17], damaging the efficiency of evolution. Redundancy is measured by
analyzing the proportion of effective mutations, and locality studies how well
genotypic neighbors correspond to phenotypic neighbors. Standard GE showed
performance similar to random search [4], which motivated researchers to
propose different initialization methods [18], representations [5, 15],
genetic operators, but also to investigate grammar design [10].

Position Independent Grammatical Evolution ($\pi$GE) [19] uses a different
representation and mapping mechanism that removes the positional dependency
that exists in GE. Each codon of the genotype contains two values, nont that
consists of the non-terminal and rule that states the rule index to expand.
This method improves performance compared with GE [20]. However, another study
found that this method also suffers from poor locality [7].

SGE [8, 5] proposed a new representation for the genotype and variation
operators, which resulted in better performance and fewer issues regarding
locality and redundancy [6, 7]. The genotype comprises a list of integers (one
per grammar non-terminal), with each integer corresponding to the index of a
production rule. This structure allows mutation to occur inside the same non-
terminal and crossover to exchange the list of derivation options. Another
advantage of this proposal is that only valid solutions are allowed. SGE
imposes a depth limit on solutions. Once the limit is surpassed, only non-
recursive productions are chosen, forcing individuals to consolidate into a
valid genotype.

Probabilistic Structured Grammatical Evolution (PSGE) [15] is a recent
proposal to SGE that uses a probabilistic grammar, namely a Probabilistic
Context-Free Grammar (PCFG), to bias the search, and where codons of the
genotype are floats. Each grammar production rule has a probability of being
selected. These probabilities change based on the frequency of expansion of
that rule on the best individual. If the rule is not expanded, the probability
decreases. This proposal performed better or similarly compared to SGE and
outperformed GE in all problems. The evolved grammar also provides information
about the features more relevant to the problem [21]. Another probabilistic
grammar approach is Co-evolutionary Probabilistic Structured Grammatical
Evolution (Co-PSGE) [22]. In this method, each individual has a PCFG, which
may suffer mutation to the probabilities values. This approach also showed
similar or better performance than SGE.

###  2.2 Adaptive mutation rate

Although most works in the literature use a static parameter for mutation and
crossover rates, research shows that dynamic parameters may improve the search
and introduce more diversity to the population. Adaptive mutations have been
widely proposed in the literature to tackle some of the issues that EAs
present.

Self-adaptive Gaussian mutation has been widely used by Evolution Strategies
(ES) [23] and adapted into EAs [24]. During the evolutionary process, the
mutation rate varies, suffering a Gaussian mutation. This approach achieves
better results than standard mutation and performs similarly to ES. Teo [25]
studied a self-adaptive Gaussian mutation operator for the Generalized
Generation Gap (G3) algorithm, and it outperformed the standard algorithm in
two of the four problems tested.

Other approaches consider individuals’ fitness when adapting the mutation
probability. Libelli et al. [26] and Lis [27] approaches showed better
performance than classical Genetic Algorithm (GA).

Adaptive Mutation Probability Genetic Algorithm (APmGA) [28] dynamically
adjusts the mutation probability during the evolutionary process based on the
variations of the population entropy between the current and previous
generations.

Salinas et al. [29] proposed an EA where operators are GP trees. In each
generation, the probability of an operator increases or decreases based on the
individual’s performance after the operator.

Gomez proposed Hybrid Adaptive Evolutionary Algorithm (HAEA) [30] to adapt the
operator probabilities during evolution. Each individual encodes its genetic
rates. The probabilities by a random value can change according to the fitness
of the offspring. This algorithm inspired other proposals that showed that
although for some problems there are no significant improvements in
performance, the algorithm can obtain similar results without the need to pre-
tuning, needing less computational time [31, 32, 33].

Coelho et al. [34] presented a new hybrid self-adaptive algorithm based on ES
guided by neighborhood structures and tested for combinatorial contains
mutation probabilities, and the second contains integer values that control
the strength of the disturbance. The results were similar to the other
approaches tested and showed that the adaptive mutation could escape local
optima and balance exploration and exploitation.

To our knowledge, there is only one adaptive mutation rate parameter proposal
in GE. Fagan et al. [35] propose Fitness Reactive Mutation (FRM), an adaptive
mutation that increases the mutation rate in case a fitness plateau is reached
to diversify the population and decreases when a new optimum is found, using
increments/decrements of 0.01. The approach found similar results as the
fixing mutation rate.

###  2.3 Grammar-design

It is possible to design grammar to produce syntactically constrained
solutions or to incorporate domain knowledge by biasing the grammar. The
grammar’s design can significantly impact the search of GE [11, 12, 13].

Miguel Nicolau [11] proposes a method to reduce the number of non-terminals of
the grammar. The authors compare standard GE with a standard and a reduced
grammar and showed an empirical increase in performance. This work motivated a
study with different types of grammars. This study shows that recursion-
balanced grammar could also improve performance [10].

A recent work by Dick et al. [12] showed that GE is more sensitive to grammar
design than Context-Free Grammar Genetic Programming (CFG-GP). The results
suggest that CFG-GP is more sensitive to parameter tuning than grammar design.

Hemberg et al. [13] compared GEs with a depth-first mapping mechanism that
uses three grammars: infix (standard GE), prefix and postfix. The results
showed that different grammars can improve performance, although the authors
report no significant differences.

Grammatical Evolution by Grammatical Evolution ((GE)2) uses two grammars, the
universal and the solution grammar. The universal grammar describes the rules
to construct the solution grammar. The rules are used to map the individuals
and can evolve towards biasing the search space. Results showed that the
evolved grammars presented some bias towards some non-terminal symbols.

Manzoni et al. [36] showed theoretically that different grammars of equal
quality impact the performance of (1+1)-EA. The structure of the grammar is
problem dependent but can favor the search. A mutation operator that modifies
the probability of selecting the grammar rules was also proposed.

##  3 Adaptive Facilitated Mutation

In biology, organisms have evolved to canalize the rate and effect of
mutations on the phenotype. Gerhart et al. [9] propose that organisms adapt to
new environments through regulatory changes that enable or disable pre-
existing conserved components. This variation increases the probability of
viable genetic mutation since core components remain unaffected by these
regulatory changes. In sum, the modularity, adaptability, and compartmentation
of genetic material in organisms allow facilitated variation through
regulatory change.

Facilitated Mutation [37] (FM) is a biologically inspired mutation mechanism
that aims to replicate the benefits of facilitated variation. This mechanism
leverages the grammar’s inherent compartmentation to regulate the mutation’s
frequency and destructiveness. With this mutation, each non-terminal has a
different mutation probability, rather than the single mutation probability
traditionally used by grammar-based EA.

In this work, we propose Adaptive Facilitated Mutation (AFM), an extension of
FM using an adaptive mutation array that removes the need to set a mutation
probability for each non-terminal manually. Each individual carries a mutation
array containing these mutation probabilities, illustrated in Figure 1. For
each individual in the initial population, the mutation array is initialized
using a specified starting mutation probability.

Figure 1: In SGE, all genes are mutated based on a single probability
(pictured left). Adaptive Facilitated Mutation uses a mutation array for each
individual containing a probability for each non-terminal (pictured right).

In each generation, all individual’s probabilities are adjusted using a random
value sampled from a Gaussian distribution with $N(0,\sigma)$, where $\sigma$
is a configurable parameter. Figure 2 illustrates the evolution of the
mutation array.

Figure 2: Example of the first perturbation to the mutation probability array
of an individual. This array is subsequently updated using a value sampled
from a Gaussian distribution. This distribution is centered at 0 using a
configurable standard deviation $\sigma$.

AFM only complements mutation operators by refining the frequency and impact
of mutation. Once the mechanism determines which non-terminals to mutate,
other operators should be used to alter the genotype within the defined scope.
During crossover, the offspring individual inherits the mutation array from
its fittest parent. More sophisticated inheritance mechanisms may improve this
approach further, but we opted for a simple strategy to validate the approach.
In this work, we use AFM for SGE and its variants, but this method is
compatible with any grammar-based GP algorithms where it is possible to tie
each codon to a corresponding non-terminal.

###  3.1 Grammar Design For Adaptive Facilitated Mutation

Since AFM leverages grammar structure, a purposefully designed grammar may
enhance the method’s performance. We hypothesize that AFM is more effective in
grammars with multiple non-terminals containing related symbols. Additionally,
a larger number of non-terminals may improve performance by enabling finer
turning of mutation probabilities through a more detailed mutation array. Non-
terminals commonly group symbols based on semantic similarity. For example, a
grammar may use a non-terminal for all operators with a single expansion
combining operators related to trigonometry (i.e., $sin$, $cos$) and the power
function (i.e., $square$, $sqrt$). These same symbols can be grouped into
several non-terminals based on function rather than semantics. Following this
reasoning, trigonometric operations would be grouped in a specific non-
terminal and the power function in a separate one. In Figure 3, we illustrate
how a grammar can be extended into a Function Grouped grammar.

$\displaystyle{<}\text{start}{>}::=$ $\displaystyle{<}\text{expr}{>}$
$\displaystyle{<}\text{expr}{>}::=$
$\displaystyle{<}\text{expr}{>}{<}\text{op}{>}{<}\text{expr}{>}\,|$
$\displaystyle{<}\text{pre\\_op}{>}{(}{<}\text{expr}{>}{)}\,|$
$\displaystyle{<}\text{var}{>}$ $\displaystyle{<}\text{op}{>}::=$
$\displaystyle{+}|{-}|{*}|{/}$ $\displaystyle{<}\text{pre\\_op}{>}::=$
$\displaystyle{sin}|{cos}|{sqrt}|{square}$ $\displaystyle{<}\text{var}{>}::=$
$\displaystyle{1.0}|{x[n]}$ (a) Initial Grammar.

$\displaystyle{<}\text{start}{>}::=$ $\displaystyle{<}\text{expr\\_var}{>}$
$\displaystyle{<}\text{expr\\_var}{>}::=$
$\displaystyle{<}\text{expr}{>}|{<}\text{var}{>}$
$\displaystyle{<}\text{expr}{>}::=$
$\displaystyle{<}\text{expr\\_var}{>}{<}\text{op}{>}{<}\text{expr\\_var}{>}\,|$
$\displaystyle{<}\text{pre\\_op}{>}{(}{<}\text{expr\\_var}{>}{)}$
$\displaystyle{<}\text{op}{>}::=$ $\displaystyle{+}|{-}|{*}|{/}$
$\displaystyle{<}\text{pre\\_op}{>}::=$
$\displaystyle{<}\text{trig\\_op}{>}|{<}\text{pow\\_op}{>}$
$\displaystyle{<}\text{trig\\_op}{>}::=$ $\displaystyle{sin}|{cos}$
$\displaystyle{<}\text{pow\\_op}{>}::=$ $\displaystyle{sqrt}|{square}$
$\displaystyle{<}\text{var}{>}::=$ $\displaystyle{1.0}|{x[n]}$ (b) Function
Grouped Grammar.

Figure 3: A grammar of semantically grouped non-terminals (3(a)) can be re-
constructed based on functional groups (3(b)). This procedure extends the
grammar and possibly improves performance when using Adaptive Facilitated
Mutation.

##  4 Experimental Setup

We use PSGE [15] in all experiments as it has equal or superior performance to
SGE in the selected tasks. We compare the algorithm using Standard Mutation
(SM) and Facilitated Mutation (FM). Additionally, we compare the Standard
Grammar (SG) and the Function Grouped Grammar (FG) that follows the principles
outlined in Section 3 (complete grammar shown in Figure 4).

We evaluate the performance of our method in three popular symbolic regression
GP benchmarks, the Quartic polynomial, the Pagie polynomial, and Boston
Housing [38, 39]. The Quartic polynomial is defined by the mathematical
expression shown in Equation 1. The function is sampled in the interval [-1,
1] with a step of 0.1.

$x[0]^{4}+x[0]^{3}+x[0]^{2}+x[0]$ (1)

The Pagie polynomial is known to be a more difficult symbolic regression
benchmark (Equation 2. The outputs are computed in the interval: $-5\leq
x[0],x[1]\leq 5.4$, with step size $0.4$. In Pagie, both features fall in the
same range of values.

$\frac{1}{1+x[0]^{-4}}+\frac{1}{1+x[1]^{-4}}$ (2)

$\displaystyle{<}\text{start}{>}::=$ $\displaystyle{<}\text{expr\\_var}{>}$
$\displaystyle{<}\text{expr\\_var}{>}::=$
$\displaystyle{<}\text{expr}{>}|{<}\text{var}{>}$
$\displaystyle{<}\text{expr}{>}::=$ $\displaystyle{<}\text{expr\\_op}{>}|$
$\displaystyle{<}\text{pre\\_op}{>}{(}{<}\text{expr\\_var}{>}{)}$
$\displaystyle{<}\text{expr\\_op}{>}::=$
$\displaystyle{<}\text{expr}{>}{<}\text{op}{>}{<}\text{expr}{>}\,|$
$\displaystyle{(}{<}\text{expr}{>}{<}\text{op}{>}{<}\text{expr}{>}{)}\,$
$\displaystyle{<}\text{op}{>}::=$ $\displaystyle{+}|{-}|{*}|{/}$
$\displaystyle{<}\text{pre\\_op}{>}::=$ $\displaystyle{<}\text{trig\\_op}{>}$
$\displaystyle|{<}\text{exp\\_log\\_op}{>}|{inv}$
$\displaystyle{<}\text{trig\\_op}{>}::=$ $\displaystyle{sin}|{cos}$
$\displaystyle{<}\text{exp\\_log\\_op}{>}::=$ $\displaystyle{exp}|{log}$
$\displaystyle{<}\text{var}{>}::=$ $\displaystyle{1.0}|{x[n]}$ Figure 4:
Function Grouped Grammar used in the experiments. Number $x[n]$ terminals
change to match the problem: 1 for Quartic, 2 for Pagie, 12 for Boston Housing

The third benchmark is Boston Housing [40]. This is a predictive modeling
problem, where one needs to build a model to predict the price of Boston
houses based on 13 features. There are 506 instances split into 90% for
training and 10% for test. The 13 features that compose the dataset are
heterogeneous regarding their intervals, with ranges varying from $0\leq
x[3]\leq 1$ to $0.32\leq x[11]\leq 396.9$. This diversity likely reduces the
effectiveness of our approach as all features are grouped in the $<var>$ non-
terminal, making it difficult for the algorithm to distinguish between
different types of features with different mutation rates. It is possible to
adjust the grammar using expert knowledge, separating the variables into non-
terminals based on orders of magnitude or function. This type of grammar
design, while possibly effective, is outside the scope of this work.

The fitness functions used to evaluate the individuals consider the
minimization of the Root Relative Squared Error (RRSE) between the
individual’s solution and the target on a data set.

Table 1 summarizes the parameters used in the experiments. All problems use
the same population size, mutation and crossover rates, tournament size, max
depth, and number of generations. In preliminary experimentation, we trialed
four parameters for the AFM’s Gaussian perturbation:
$\sigma=[0.001,0.0025,0.005,0.01]$ We found that facilitated mutations’
$\sigma$ parameter benefited from tuning when moved to different tasks.
Consequently, each experiment uses the best $\sigma$ value for the
corresponding problem. We repeat all experiments 100 times to investigate
statistically meaningful differences between the approaches.

Table 1: Parameters used in experiments for Quartic, Pagie, and Boston Housing. Parameters | Quartic | Pagie |  | Boston
---
Housing
Population Size | 1000
Generations | 100
Elitism | 10%
Mutation | Gaussian N(0,0.5)
Mutation Probability | 10%
Adaptive Facilitated Mutation $\sigma$ | 0.0025 | 0.001
Crossover Probability | 90%
Tournament Size | 3
Max Depth | 10

##  5 Results

This section presents the results obtained for each problem in terms of the
Mean Best Fitness (MBF) of 100 repetitions. We statistically compare the
different approaches using the Mann-Whitney test with Bonferroni correction
with a significance level $\alpha=0.05$.

Figure 5: Plot shows the mean best fitness of 100 runs for the Quartic
polynomial.

Figure 5 shows the results for the Quartic polynomial. At the end of
evolution, all approaches achieve a similar result in this problem, except for
SM+SG, which performs significantly worse than the others (see Table 2 for
statistical analysis). Quartic likely has a local optimum that evolution
cannot escape without our proposed mechanism. While SM+FG and AFM+SG
statistically outperform the baseline, AFM+FG’s advantages are only empirical.
We hypothesize that AFM+FG would join the other two methods with additional
repetitions, creating two tiers of solution quality in Quartic.

Table 2: P-value for Mann-Whitney Statistical Tests using Bonferroni correction with significance level $\alpha=0.05$ for Quartic polynomial. Bold indicates that the method in the corresponding row is statistically superior. | Quartic
---
polynomial
SM+SG | SM+FG | AFM+SG | AFM+FG
SM+SG |  |  |  |
SM+FG | 0.048 |  |  |
AFM+SG | 0.004 | 0.356 |  |
AFM+FG | 0.072 | 0.760 | 0.332 |
Figure 6: Plot shows the mean best fitness of 100 runs for the Pagie
polynomial.

Figure 6 shows the MBF for the Pagie polynomial across 100 generations.
Looking at the results, one can see that the advantages of the proposed
approaches are clear for this problem, particularly the FG grammar. From
generation 25 until termination, both AFM solutions have a better MBF than
their SM counterparts. While AFM is empirically superior to SM, FG amplifies
the benefits of this approach.

The statistical analysis results (shown in Table 3) reveal that AFM+FG
outperforms all SG approaches while SM+FG is similar. These results suggest
that our proposed grammar may be marginally superior, but considerable
benefits come from combining it with the appropriate mutation.

Table 3: P-value for Mann-Whitney Statistical Tests using Bonferroni correction with significance level $\alpha=0.05$ for Pagie polynomial. Bold indicates that the method in the corresponding row is statistically superior. | Pagie
---
polynomial
SM+SG | SM+FG | AFM+SG | AFM+FG
SM+SG |  |  |  |
SM+FG | 0.965 |  |  |
AFM+SG | 0.426 | 0.400 |  |
AFM+FG | 0.002 | 0.001 | 0.001 |

(a) Training

(b) Test

Figure 7: Plot shows the mean best fitness of 100 runs for the Boston Housing
dataset.

Finally, in Figure 7 we show the MBF for the Boston Housing Training and Test.
A brief perusal of the results indicates that AFM+FG achieves the best MBF in
Test. When comparing the mutations, AFM generalizes better as it maintains a
similar performance between train and test.

Note that SM+SG appears to be worst at generalizing to the test data. This is
most evident when comparing SM+SG with both AFM approaches, as the differences
between these methods are noticeably larger when moved to test data. Despite
the highlighted differences, statistical tests for Boston Housing Training and
Test reveal no statistical differences (Tables 4 and 5). Given that FG cannot
account for the feature variety of Boston Housing, it is remarkable that the
approach still achieves competitive results, especially in the test data. It
is possible that a specifically designed FG that uses expert knowledge to
group the features based on function could achieve even better and
statistically significant results.

Table 4: P-value for Mann-Whitney Statistical Tests using Bonferroni correction with significance level $\alpha=0.05$ for Boston Housing Training. Bold indicates that the method in the corresponding row is statistically superior. | Boston Housing
---
Training
SM+SG | SM+FG | AFM+SG | AFM+FG
SM+SG |  |  |  |
SM+FG | 0.530 |  |  |
AFM+SG | 0.263 | 0.763 |  |
AFM+FG | 0.531 | 0.361 | 0.548 |
Table 5: P-value for Mann-Whitney Statistical Tests using Bonferroni correction with significance level $\alpha=0.05$ for Boston Housing Test. Bold indicates that the method in the corresponding row is statistically superior. | Boston Housing
---
Test
SM+SG | SM+FG | AFM+SG | AFM+FG
SM+SG |  |  |  |
SM+FG | 0.929 |  |  |
AFM+SG | 0.377 | 0.333 |  |
AFM+FG | 0.364 | 0.403 | 0.963 |

##  6 Conclusion

In this paper, we propose AFM, a mutation method that leverages grammar-based
GP’s properties to replicate natural evolutionary phenomena. This approach
divides the single mutation probability commonly found in such approaches into
a mutation array, where each grammar non-terminal has a corresponding mutation
rate. Each individual has a mutation array that co-evolves with the genetic
code. A randomly sampled value from a Gaussian distribution adjusts the
mutation rates of all individuals in each generation.

We also propose a grammar-design approach, Function Grouped Grammars, to
enhance the effectiveness of the mutation proposed. Function Grouped Grammars
organize non-terminals based on functional similarity rather than the semantic
similarity common in the field. We compare our proposals with a baseline
(standard mutation and standard grammar) and find that, when combined, our
approaches are statistically superior or similar to the baseline in three
relevant GP benchmarks.

This approach still requires parameter tuning, one of the problems tackled by
the literature by proposing adaptive mutations. However, few of these
approaches consider different values for different symbols [34]. This work
shows that an adaptive mutation rate can be beneficial for search and grammar
design amplifies these benefits.

###  6.1 Future Work

The results of our experiments are promising, but additional tests are
essential to validate our approach further. The benchmarks addressed are
relevant in GP, but a more extensive (and varied) set of benchmarks could
bring meaningful insights into the general applicability of the FG+FM.

Regarding AFM, we use a fixed starting mutation rate in all experiments. Our
method leverages adaptability as a tool for improved evolution, but we did not
investigate the potential of AFM as a replacement for mutation rate tuning.
Further development of AFM may also lead to an implementation that does not
rely on Gaussian distributions and $sigma$ tuning for perturbations. In the
future, we want to explore alternatives where AFM is a competitive, parameter-
less alternative to standard mutation. Another line of work is to experiment
with different inheritance mechanisms during crossover. In this work, the
offspring inherited the array of probabilities from the most fitted parent. It
would be interesting to explore the random selection of the parent that passes
the array or the application of the existing SGE crossover to the mutation
arrays of the parents. Such approaches would better preserve the advantages of
co-evolution, possibly improving results.

Function Grouped Grammars can also be further investigated. In this work, we
still apply this idea conservatively. Considering the FG used in experiments,
it is still possible to separate constants from variables and commutative
operators from non-commutative operators. Barring small details, we use the
same grammar for all tasks. Applying the same approach to problems requiring
more complex grammar would be interesting, as such problems yield more
opportunities for function grouping.

Finally, the applicability of these ideas to different, compatible grammar-
based GP must also be investigated. While this works focuses on PSGE, the same
approach is easily applicable in SGE [5], Co-PSGE [22], and even more
different systems like $\pi$GE [19]. Any grammar-based approach where genes
are tied to a non-terminal may benefit from FG+FM.

## Acknowledgments

This work was funded by FEDER funds through the Operational Programme
Competitiveness Factors - COMPETE and national funds by FCT - Foundation for
Science and Technology (POCI-01-0145-FEDER-029297, CISUC - UID/CEC/00326/2020)
and within the scope of the project A4A: Audiology for All
(CENTRO-01-0247-FEDER-047083) financed by the Operational Program for
Competitiveness and Internationalisation of PORTUGAL 2020 through the European
Regional Development Fund. The first author is funded by FCT, Portugal, under
the grant UI/BD/151053/2021 and the second under the grant 2022.10174.BD.

## References

  * [1] Whigham, P.A., Science, D.O.C.: Grammatically-based genetic programming (1995)
  * [2] Ryan, C., Collins, J., Neill, M.O.: Grammatical evolution: Evolving programs for an arbitrary language. In: Lecture Notes in Computer Science, pp. 83–96. Springer Berlin Heidelberg (1998)
  * [3] McKay, R.I., Hoai, N.X., Whigham, P.A., Shan, Y., O’Neill, M.: Grammar-based genetic programming: a survey. Genetic Programming and Evolvable Machines 11(3-4), 365–396 (May 2010)
  * [4] Whigham, P.A., Dick, G., Maclaurin, J., Owen, C.A.: Examining the ”best of both worlds” of grammatical evolution. In: Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation. ACM (Jul 2015)
  * [5] Lourenço, N., Assunção, F., Pereira, F.B., Costa, E., Machado, P.: Structured grammatical evolution: A dynamic approach. In: Handbook of Grammatical Evolution, pp. 137–161. Springer International Publishing (2018)
  * [6] Lourenço, N., Ferrer, J., Pereira, F.B., Costa, E.: A comparative study of different grammar-based genetic programming approaches. In: Lecture Notes in Computer Science, pp. 311–325. Springer International Publishing (2017)
  * [7] Medvet, E.: A comparative analysis of dynamic locality and redundancy in grammatical evolution. In: Lecture Notes in Computer Science, pp. 326–342. Springer International Publishing (2017)
  * [8] Lourenço, N., Pereira, F.B., Costa, E.: Unveiling the properties of structured grammatical evolution. Genetic Programming and Evolvable Machines 17(3), 251–289 (Feb 2016)
  * [9] Gerhart, J., Kirschner, M.: The theory of facilitated variation. Proceedings of the National Academy of Sciences 104(suppl_1), 8582–8589 (2007), https://www.pnas.org/doi/abs/10.1073/pnas.0701035104
  * [10] Nicolau, M., Agapitos, A.: Understanding grammatical evolution: Grammar design. In: Handbook of Grammatical Evolution, pp. 23–53. Springer International Publishing (2018)
  * [11] Nicolau, M.: Automatic grammar complexity reduction in grammatical evolution. In: The 3rd Grammatical Evolution Workshop: A workshop of the 2004 Genetic and Evolutionary Computation Conference (GECCO-2004), Seattle, Washington, USA, 26-30 June 2004. Seattle, Washington, USA (26-30 2004)
  * [12] Dick, G., Whigham, P.A.: Initialisation and grammar design in grammar-guided evolutionary computation. In: Proceedings of the Genetic and Evolutionary Computation Conference Companion. p. 534–537. GECCO ’22, Association for Computing Machinery, New York, NY, USA (2022), https://doi.org/10.1145/3520304.3529051
  * [13] Hemberg, E.: Pre-, in-and postfix grammars for symbolic regression in grammatical evolution (2008)
  * [14] Ryan, C., O'Neill, M., Collins, J. (eds.): Handbook of Grammatical Evolution. Springer International Publishing (2018)
  * [15] Megane, J., Lourenco, N., Machado, P.: Probabilistic structured grammatical evolution. In: 2022 IEEE Congress on Evolutionary Computation (CEC). IEEE (Jul 2022), https://doi.org/10.1109/cec55065.2022.9870397
  * [16] Thorhauer, A.: On the non-uniform redundancy in grammatical evolution. In: Parallel Problem Solving from Nature – PPSN XIV, pp. 292–302. Springer International Publishing (2016)
  * [17] Rothlauf, F., Oetzel, M.: On the locality of grammatical evolution. In: Lecture Notes in Computer Science, pp. 320–330. Springer Berlin Heidelberg (2006)
  * [18] Nicolau, M.: Understanding grammatical evolution: initialisation. Genetic Programming and Evolvable Machines 18(4), 467–507 (Jul 2017)
  * [19] O’Neill, M., Brabazon, A., Nicolau, M., McGarraghy, S., Keenan, P.: $\pi$grammatical evolution. In: Genetic and Evolutionary Computation – GECCO 2004, pp. 617–629. Springer Berlin Heidelberg (2004)
  * [20] Fagan, D., O’Neill, M., Galván-López, E., Brabazon, A., McGarraghy, S.: An analysis of genotype-phenotype maps in grammatical evolution. In: Lecture Notes in Computer Science, pp. 62–73. Springer Berlin Heidelberg (2010)
  * [21] Mégane, J., Lourenço, N., Machado, P.: Probabilistic grammatical evolution. In: Hu, T., Lourenço, N., Medvet, E. (eds.) Genetic Programming. pp. 198–213. Springer International Publishing, Cham (2021)
  * [22] Mégane, J., Lourenço, N., Machado, P.: Co-evolutionary probabilistic structured grammatical evolution. In: Proceedings of the Genetic and Evolutionary Computation Conference. ACM (Jul 2022), https://doi.org/10.1145/3512290.3528833
  * [23] Beyer, H., Schwefel, H.: Evolution strategies – a comprehensive introduction. Natural Computing 1, 3–52 (2004)
  * [24] Hinterding, R.: Gaussian mutation and self-adaption for numeric genetic algorithms. In: Proceedings of 1995 IEEE International Conference on Evolutionary Computation. IEEE
  * [25] Teo, J.: Self-adaptive mutation for enhancing evolutionary search in real-coded genetic algorithms. In: 2006 International Conference on Computing & Informatics. IEEE (Jun 2006)
  * [26] Libelli, S.M., Alba, P.: Adaptive mutation in genetic algorithms. Soft Computing 4(2), 76–80 (Jul 2000)
  * [27] Lis, J.: Genetic algorithm with the dynamic probability of mutation in the classification problem. Pattern Recognition Letters 16(12), 1311–1320 (Dec 1995)
  * [28] Stark, N., Minetti, G.F., Salto, C.: A new strategy for adapting the mutation probability in genetic algorithms (2012)
  * [29] Cruz-Salinas, A.F., Perdomo, J.G.: Self-adaptation of genetic operators through genetic programming techniques. In: Proceedings of the Genetic and Evolutionary Computation Conference. ACM (jul 2017)
  * [30] Gómez, J.: Self adaptation of operator rates in evolutionary algorithms. In: GECCO (2004)
  * [31] Gómez, J., León, E.: On the class of hybrid adaptive evolutionary algorithms (chavela). Natural Computing 20(3), 377–394 (Feb 2021)
  * [32] Montero, E., Riff, M.C.: Calibrating strategies for evolutionary algorithms. In: 2007 IEEE Congress on Evolutionary Computation. IEEE (Sep 2007)
  * [33] Montero, E., Riff, M.C.: Self-calibrating strategies for evolutionary approaches that solve constrained combinatorial problems. In: Lecture Notes in Computer Science, pp. 262–267. Springer Berlin Heidelberg
  * [34] Coelho, V.N., Coelho, I.M., Souza, M.J.F., Oliveira, T.A., Cota, L.P., Haddad, M.N., Mladenovic, N., Silva, R.C.P., Guimarães, F.G.: Hybrid self-adaptive evolution strategies guided by neighborhood structures for combinatorial optimization problems. Evolutionary Computation 24(4), 637–666 (Dec 2016)
  * [35] Fagan, D., Hemberg, E., Nicolau, M., O'Neill, M., McGarraghy, S.: Towards adaptive mutation in grammatical evolution. In: Proceedings of the fourteenth international conference on Genetic and evolutionary computation conference companion - GECCO Companion '12. ACM Press (2012)
  * [36] Manzoni, L., Bartoli, A., Castelli, M., Goncalves, I., Medvet, E.: Specializing context-free grammars with a (1 + 1)-EA. IEEE Transactions on Evolutionary Computation 24(5), 960–973 (oct 2020)
  * [37] Tiso, S., Carvalho, P., Lourenço, N., Machado, P.: Structured mutation inspired by evolutionary theory enriches population performance and diversity (2023), https://arxiv.org/abs/2302.00559
  * [38] White, D.R., McDermott, J., Castelli, M., Manzoni, L., Goldman, B.W., Kronberger, G., Jaśkowski, W., O’Reilly, U., Luke, S.: Better GP benchmarks: community survey results and proposals. Genetic Programming and Evolvable Machines 14(1), 3–29 (Dec 2012)
  * [39] McDermott, J., Jong, K.D., O'Reilly, U., White, D.R., Luke, S., Manzoni, L., Castelli, M., Vanneschi, L., Jaskowski, W., Krawiec, K., Harper, R.: Genetic programming needs better benchmarks. In: Proceedings of the fourteenth international conference on Genetic and evolutionary computation conference - GECCO '12. ACM Press (2012)
  * [40] Harrison Jr, D., Rubinfeld, D.L.: Hedonic housing prices and the demand for clean air. Journal of environmental economics and management 5(1), 81–102 (1978)


# On saturation of the discrepancy principle for nonlinear Tikhonov
regularization in Hilbert spaces

Qinian Jin  Mathematical Sciences Institute, Australian National University,
Canberra, ACT 2601, Australia<EMAIL_ADDRESS>
###### Abstract.

In this paper we revisit the discrepancy principle for Tikhonov regularization
of nonlinear ill-posed problems in Hilbert spaces and provide some new and
improved saturation results under less restrictive conditions, comparing with
the existing results in the literature.

###### Key words and phrases:

Nonlinear ill-posed problems, Tikhonov regularization, the discrepancy
principle, saturation

##  1\. Introduction

Let $X$ and $Y$ be two Hilbert spaces whose inner products and norms are
denoted by the same notations $\langle\cdot,\cdot\rangle$ and $\|\cdot\|$. Let
$F:\mbox{dom}(F)\subset X\to Y$ be a nonlinear operator with domain
$\mbox{dom}(F)$. We consider the equation

$F(x)=y$ (1)

and assume (1) is solvable in $\mbox{dom}(F)$, i.e. $y\in\mbox{Ran}(F)$, the
range of $F$. Since the solution of (1) may not be unique, it is necessary to
select a right one. To this end, we make a first guess $x^{*}\in\mbox{dom}(F)$
and try to find a solution closest to $x^{*}$. We say
$x^{\dagger}\in\mbox{dom}(F)$ is a $x^{*}$-minimum-norm solution of (1) if
$F(x^{\dagger})=y$ and

$\|x^{\dagger}-x^{*}\|=\min_{x\in\mbox{dom}(F)}\\{\|x-x^{*}\|:F(x)=y\\}.$

We consider the situation that (1) is ill-posed in the sense that its solution
does not depend continuously on the data. The exact data $y$ is usually
unavailable in practical applications. Instead we only have a noisy data
$y^{\delta}$ obtained by measurement which satisfies

$\|y^{\delta}-y\|\leq\delta$ (2)

with a known noise level $\delta>0$. One of the important questions is then to
reconstruct $x^{\dagger}$ by using $y^{\delta}$. This requires the
regularization techniques. The most famous method is Tikhonov regularization
in which the solution of the minimization problem

$\min_{x\in\mbox{dom}(F)}\left\\{\|F(x)-y^{\delta}\|^{2}+\alpha\|x-x^{*}\|^{2}\right\\}$
(3)

is used to approximate $x^{\dagger}$, where $\alpha>0$ is the so-called
regularization parameter. The existence of $x^{*}$-minimum-norm solution and
the regularized solution $x_{\alpha}^{\delta}$ can be guaranteed when $F$ is
weakly closed; see [3]. The solution of (3) may not be unique in general, we
will use $x_{\alpha}^{\delta}$ to denote any solution of (3). When
$x_{\alpha}^{\delta}$ is used to approximate $x^{\dagger}$, the accuracy
depends crucially on the choice of the regularization parameter $\alpha$ which
can be selected by either a priori or a posteriori rules ([3, 8, 11, 12]). One
of the most prominent a posteriori parameter choice rule is Morozov’s
discrepancy principle which can be stated as follows.

######  Rule 1 (The discrepancy principle).

Let $\tau\geq 1$ be a given number. Define $\alpha(\delta,y^{\delta})>0$ to be
the number such that
$\|F(x_{\alpha(\delta,y^{\delta})}^{\delta})-y^{\delta}\|=\tau\delta$.

One may refer to [7] for a discussion on the existence of
$\alpha(\delta,y^{\delta})$ satisfying Rule 1. With
$\alpha(\delta,y^{\delta})$ chosen by Rule 1, it has been shown in [3] that
$x_{\alpha(\delta,y^{\delta})}^{\delta}\to x^{\dagger}$ as $\delta\to 0$ and
if $F$ is Fréchet differentiable in a neighborhood of $x^{\dagger}$ whose
Fréchet derivative $F^{\prime}(x)$ is Lipschitz continuous with constant $L$
and if the source condition $x^{\dagger}-x^{*}=F^{\prime}(x^{\dagger})^{*}u$
holds with $L\|u\|<1$, then

$\displaystyle\sup\left\\{\|x_{\alpha(\delta,y^{\delta})}^{\delta}-x^{\dagger}\|:\|y^{\delta}-y\|\leq\delta\right\\}=O(\delta^{1/2})$
(4)

for the worst case error. It is natural to ask if the rate in (4) can be
improved under a higher order source condition
$x^{\dagger}-x^{*}=\mbox{Ran}((F^{\prime}(x^{\dagger})^{*}F^{\prime}(x^{\dagger}))^{\nu})$
with $\nu>1/2$. For compact linear operator $F$, Groetsch showed in [4, 5]
that if the rate $O(\delta^{1/2})$ in (4) is replaced by $o(\delta^{1/2})$
then it must have $x^{*}=x^{\dagger}$. Thus, for a general choice of the
initial guess $x^{*}$, the rate in (4) is the best possible in the worst case
scenario which is known as the saturation of the discrepancy principle for
linear Tikhonov regularization in Hilbert spaces. In [10] Scherzer gave a
saturation result on Rule 1 for nonlinear Tikhonov regularization in Hilbert
spaces. His proof relies on several conditions on the nonlinear operator $F$
some of which seem restrictive. Furthermore, the argument in [10] requires
$\tau$ to be sufficiently large which seems unnatural. In this note we will
provide new and improved saturation results for Rule 1 for nonlinear Tikhonov
regularization (3) in Hilbert spaces under less restrictive conditions.

##  2\. Main results

The saturation result proved in [10] on Rule 1 for nonlinear Tikhonov
regularization in Hilbert spaces is based on the following technical
conditions on the operator $F$, where $B_{\rho}(x^{\dagger}):=\\{x\in
X:\|x-x^{\dagger}\|<\rho\\}$.

######  Assumption 1.

There is $\rho>0$ such that $B_{\rho}(x^{\dagger})\subset\emph{dom}(F)$ and
$F$ is Fréchet dfferentiable on $B_{\rho}(x^{\dagger})$. Moreover, there
exists $\kappa_{0}\geq 0$ such that for any $x,z\in B_{\rho}(x^{\dagger})$
there is a bounded linear operator $R(x,z):X\to X$ such that

$F^{\prime}(x)=F^{\prime}(z)R(x,z)\quad\mbox{and}\quad\|I-R(x,z)\|\leq\kappa_{0}\|x-z\|.$

######  Assumption 2.

There exists $\kappa_{1}\geq 0$ such that for any $x,z\in
B_{\rho}(x^{\dagger})$ there is a bounded linear operator $S(x,z):Y\to Y$ such
that

$F^{\prime}(x)=S(x,z)F^{\prime}(z)\quad\mbox{and}\quad\|I-S(x,z)\|\leq\kappa_{1}\|x-z\|.$

######  Theorem 1 ([10]).

Let Assumption 1 and Assumption 2 hold. Assume that $F^{\prime}(x^{\dagger})$
is compact with infinite rank and
$\max\\{\kappa_{0},\kappa_{1}\\}\|x^{*}-x^{\dagger}\|$ is sufficiently small.
Consider Tikhonov regularization (3) and let $\alpha(\delta,y^{\delta})$ be
determined by Rule 1 with sufficiently large $\tau>1$. If

$\sup\left\\{\|x_{\alpha(\delta,y^{\delta})}^{\delta}-x^{\dagger}\|:\|y^{\delta}-y\|\leq\delta\right\\}=o(\delta^{1/2})$

as $\delta\to 0$, then $x^{\dagger}=x^{*}$.

The proof idea of Theorem 1 given in [10] is to connect Rule 1 for the
nonlinear Tikhonov regularization (3) with a discrepancy principle for its
linearization at the sought solution $x^{\dagger}$ and then use the saturation
result of Groetsch ([4, 5]) for linear Tikhonov regularization. More
precisely, replacing $F(x)$ in (3) by its linearization $y+A(x-x^{\dagger})$
at $x^{\dagger}$, where $A:=F^{\prime}(x^{\dagger})$, one may consider the
linearized minimization problem

$\displaystyle\min_{x\in
X}\left\\{\|Ax-q^{\delta}\|^{2}+\alpha\|x-x^{*}\|^{2}\right\\},\quad
q^{\delta}:=Ax^{\dagger}+y^{\delta}-y$ (5)

whose unique minimizer is given by

$\hat{x}_{\alpha}^{\delta}=x^{\dagger}+(\alpha
I+A^{*}A)^{-1}\left(\alpha(x^{*}-x^{\dagger})+A^{*}(y^{\delta}-y)\right).$ (6)

It is then shown that if $\alpha$ is chosen by Rule 1, there must exists a
constant $c>1$ such that

$\displaystyle\|A\hat{x}_{\alpha}^{\delta}-q^{\delta}\|=\hat{\tau}\delta$ (7)

for some number $\hat{\tau}\in[1,c]$. Once this is achieved, the saturation
result on the discrepancy principle for linear Tikhonov regularization can be
used directly to conclude $x^{*}-x^{\dagger}=0$. In order to achieve the above
connection between Rule 1 for nonlinear Tikhonov regularization and a
discrepancy principle for its linearization,one needs to compare
$\|F(x_{\alpha}^{\delta})-y^{\delta}\|$ with
$\|A\hat{x}_{\alpha}^{\delta}-q^{\delta}\|$ and guarantee the existence of a
constant $c$ such that (7) holds for some $\hat{\tau}\in[1,c]$. These require
the number $\tau$ in Rule 1 to be sufficiently large and the nonlinear
operator $F$ to satisfy Assumption 1 and Assumption 2 with sufficiently small
$\max\\{\kappa_{0},\kappa_{1}\\}\|x^{*}-x^{\dagger}\|$.

In the following we will present some new and improved saturation results on
Rule 1 for nonlinear Tikhonov regularization (3) in Hilbert spaces under less
restrictive conditions. Our first result shows that the conclusion of Theorem
1 holds for Rule 1 with any $\tau>0$ under merely Assumption 1 and the
smallness of $\max\\{\kappa_{0},\kappa_{1}\\}\|x^{*}-x^{\dagger}\|$ is not
required.

######  Theorem 2.

Let Assumption 1 hold. Assume that, for the operator
$A:=F^{\prime}(x^{\dagger})$, there exists a sequence of positive numbers
$\\{\lambda_{k}\\}$ in the spectrum $\sigma(AA^{*})$ of $AA^{*}$ such that
$\lim_{k\to\infty}\lambda_{k}=0$. Consider the nonlinear Tikhonov
regularization (3) and assume that there is a constant $\tau>0$ such that, for
any small $\delta>0$ and any $y^{\delta}$ satisfying
$\|y^{\delta}-y\|\leq\delta$, a number $\alpha(\delta,y^{\delta})>0$ can be
chosen such that
$\|F(x_{\alpha(\delta,y^{\delta})}^{\delta})-y^{\delta}\|\leq\tau\delta$. If

$\displaystyle\sup\left\\{\|x_{\alpha(\delta,y^{\delta})}^{\delta}-x^{\dagger}\|:\|y^{\delta}-y\|\leq\delta\right\\}=o(\delta^{1/2})$
(8)

as $\delta\to 0$, then $x^{\dagger}=x^{*}$.

###### Proof.

Let $\alpha:=\alpha(\delta,y^{\delta})$. According to (8), we have
$\|x_{\alpha}^{\delta}-x^{\dagger}\|<\rho$ for all $y^{\delta}$ satisfying
$\|y^{\delta}-y\|\leq\delta$ provided $\delta>0$ is sufficiently small. Thus
$x_{\alpha}^{\delta}$ is an interior point of $\mbox{dom}(F)$ and consequently
there holds the Euler equation

$F^{\prime}(x_{\alpha}^{\delta})^{*}(F(x_{\alpha}^{\delta})-y^{\delta})+\alpha(x_{\alpha}^{\delta}-x^{*})=0.$

Manipulating this equation we then obtain

$\displaystyle(\alpha I+A^{*}A)(x_{\alpha}^{\delta}-x^{\dagger})$
$\displaystyle=\alpha(x^{*}-x^{\dagger})+F^{\prime}(x_{\alpha}^{\delta})^{*}(y^{\delta}-F(x_{\alpha}^{\delta}))+A^{*}A(x_{\alpha}^{\delta}-x^{\dagger})$
$\displaystyle=\alpha(x^{*}-x^{\dagger})+(F^{\prime}(x_{\alpha}^{\delta})^{*}-A^{*})(y^{\delta}-F(x_{\alpha}^{\delta}))$
$\displaystyle\quad\,-A^{*}(F(x_{\alpha}^{\delta})-y^{\delta}-A(x_{\alpha}^{\delta}-x^{\dagger}))$

which gives

$\displaystyle x_{\alpha}^{\delta}-x^{\dagger}$ $\displaystyle=(\alpha
I+A^{*}A)^{-1}\left(\alpha(x^{*}-x^{\dagger})+A^{*}(y^{\delta}-y)\right)$
$\displaystyle\quad\,+(\alpha
I+A^{*}A)^{-1}(F^{\prime}(x_{\alpha}^{\delta})^{*}-A^{*})(y^{\delta}-F(x_{\alpha}^{\delta}))$
$\displaystyle\quad\,-(\alpha
I+A^{*}A)^{-1}A^{*}(F(x_{\alpha}^{\delta})-y-A(x_{\alpha}^{\delta}-x^{\dagger})).$

By using the definition of $\hat{x}_{\alpha}^{\delta}$, we thus have

$\displaystyle x_{\alpha}^{\delta}-\hat{x}_{\alpha}^{\delta}$
$\displaystyle=-(\alpha
I+A^{*}A)^{-1}A^{*}\left(F(x_{\alpha}^{\delta})-y-A(x_{\alpha}^{\delta}-x^{\dagger})\right)$
$\displaystyle\quad\,+(\alpha
I+A^{*}A)^{-1}\left(F^{\prime}(x_{\alpha}^{\delta})^{*}-A^{*}\right)\left(y^{\delta}-F(x_{\alpha}^{\delta})\right).$
(9)

Since $x_{\alpha}^{\delta}\in B_{\rho}(x^{\dagger})$, we may use Assumption 1
to obtain

$\displaystyle x_{\alpha}^{\delta}-\hat{x}_{\alpha}^{\delta}$
$\displaystyle=-\int_{0}^{1}(\alpha
I+A^{*}A)^{-1}A^{*}A\left(R(x^{\dagger}+t(x_{\alpha}^{\delta}-x^{\dagger}),x^{\dagger})-I\right)(x_{\alpha}^{\delta}-x^{\dagger})dt$
$\displaystyle\quad\,+(\alpha
I+A^{*}A)^{-1}\left(I-R(x^{\dagger},x_{\alpha}^{\delta})^{*}\right)F^{\prime}(x_{\alpha}^{\delta})^{*}(y^{\delta}-F(x_{\alpha}^{\delta})).$

and consequently

$\displaystyle\|x_{\alpha}^{\delta}-\hat{x}_{\alpha}^{\delta}\|$
$\displaystyle\leq\int_{0}^{1}t\kappa_{0}\|x_{\alpha}^{\delta}-x^{\dagger}\|^{2}dt+\frac{1}{\alpha}\|I-R(x^{\dagger},x_{\alpha}^{\delta})\|\|F^{\prime}(x_{\alpha}^{\delta})^{*}(y^{\delta}-F(x_{\alpha}^{\delta}))\|$
$\displaystyle\leq\frac{1}{2}\kappa_{0}\|x_{\alpha}^{\delta}-x^{\dagger}\|^{2}+\kappa_{0}\|x_{\alpha}^{\delta}-x^{\dagger}\|\|x_{\alpha}^{\delta}-x^{*}\|\leq
C_{0}\|x_{\alpha}^{\delta}-x^{\dagger}\|,$

where $C_{0}:=(3\rho/2+\|x^{*}-x^{\dagger}\|)\kappa_{0}$. Therefore

$\displaystyle\|\hat{x}_{\alpha}^{\delta}-x^{\dagger}\|\leq(1+C_{0})\|x_{\alpha}^{\delta}-x^{\dagger}\|.$

By virtue of (8) we then have

$\displaystyle\sup\left\\{\|\hat{x}_{\alpha}^{\delta}-x^{\dagger}\|:\|y^{\delta}-y\|\leq\delta\right\\}=o(\delta^{1/2}).$
(10)

Based on (10) and the assumption on $\sigma(AA^{*})$, we will show that for
$\delta_{k}:=\lambda_{k}$ there exists $y^{\delta_{k}}$ satisfying
$\|y^{\delta_{k}}-y\|\leq\delta_{k}$ such that

$\delta_{k}/\alpha_{k}\to 0\quad\mbox{ as }k\to\infty,$

where $\alpha_{k}:=\alpha(\delta_{k},y^{\delta_{k}})$. To see this, let
$\\{F_{\lambda}\\}$ denote the spectral family generated by the self-adjoint
operator $AA^{*}$ and for each $k$ we set

$G_{\lambda_{k}}:=F_{3\lambda_{k}/2}-F_{\lambda_{k}/2}.$

Then $G_{\lambda_{k}}$ is an orthogonal projection on $Y$ and commutes with
$(\alpha I+AA^{*})^{-1}$. We then follow an idea from [1, 9] and consider the
perturbed data defined by

$y^{\delta_{k}}:=y-\delta_{k}G_{\lambda_{k}}z_{k},$

where

$z_{k}:=\left\\{\begin{array}[]{lll}\|G_{\lambda_{k}}A(x^{\dagger}-x^{*})\|^{-1}A(x^{\dagger}-x^{*}),&\mbox{
if }G_{\lambda_{k}}A(x^{\dagger}-x^{*})\neq 0,\\\ \mbox{arbitrary with
}\|G_{\lambda_{k}}z_{k}\|=1,&\mbox{ otherwise}.\end{array}\right.$

Since $\lambda_{k}\in\sigma(AA^{*})$, such $y^{\delta_{k}}$ is well-defined.
It is clear that $\|y^{\delta_{k}}-y\|=\delta_{k}$. Let

$\hat{x}_{\alpha_{k}}:=x^{\dagger}+\alpha_{k}(\alpha_{k}I+A^{*}A)^{-1}(x^{*}-x^{\dagger}).$

From the definitions of $\hat{x}_{\alpha_{k}}^{\delta_{k}}$ and
$\hat{x}_{\alpha_{k}}$ it follows that

$\hat{x}_{\alpha_{k}}^{\delta_{k}}-\hat{x}_{\alpha_{k}}=-\delta_{k}(\alpha_{k}I+A^{*}A)^{-1}A^{*}G_{\lambda_{k}}z_{k}.$

Hence, by using $G_{\lambda_{k}}^{2}=G_{\lambda_{k}}$, the definition of
$z_{k}$ and $\|G_{\lambda_{k}}z_{k}\|=1$, we have

$\displaystyle
2\left\langle\hat{x}_{\alpha_{k}}-x^{\dagger},\hat{x}_{\alpha_{k}}^{\delta_{k}}-\hat{x}_{\alpha_{k}}\right\rangle$
$\displaystyle=-2\alpha_{k}\delta_{k}\left\langle(\alpha_{k}I+AA^{*})^{-1}A(x^{*}-x^{\dagger}),(\alpha_{k}I+AA^{*})^{-1}G_{\lambda_{k}}z_{k}\right\rangle$
$\displaystyle=2\alpha_{k}\delta_{k}\left\langle(\alpha_{k}I+AA^{*})^{-1}G_{\lambda_{k}}A(x^{\dagger}-x^{*}),(\alpha_{k}I+AA^{*})^{-1}G_{\lambda_{k}}z_{k}\right\rangle$
$\displaystyle=2\alpha_{k}\delta_{k}\|G_{\lambda_{k}}A(x^{\dagger}-x^{*})\|\|(\alpha_{k}I+AA^{*})^{-1}G_{\lambda_{k}}z_{k}\|^{2}$
$\displaystyle\geq 0$

and

$\displaystyle\|\hat{x}_{\alpha_{k}}^{\delta_{k}}-\hat{x}_{\alpha_{k}}\|^{2}$
$\displaystyle=\delta_{k}^{2}\|(AA^{*})^{1/2}(\alpha_{k}I+AA^{*})^{-1}G_{\lambda_{k}}z_{k}\|^{2}$
$\displaystyle=\delta_{k}^{2}\int_{\lambda_{k}/2}^{3\lambda_{k}/2}\frac{\lambda}{(\alpha_{k}+\lambda)^{2}}d\|F_{\lambda}G_{\lambda_{k}}z_{k}\|^{2}$
$\displaystyle\geq\delta_{k}^{2}\frac{\lambda_{k}/2}{(\alpha_{k}+3\lambda_{k}/2)^{2}}\|G_{\lambda_{k}}z_{k}\|^{2}$
$\displaystyle=\delta_{k}^{2}\frac{\lambda_{k}/2}{(\alpha_{k}+3\lambda_{k}/2)^{2}}.$

Therefore, by using the identity

$\displaystyle\|\hat{x}_{\alpha_{k}}^{\delta_{k}}-x^{\dagger}\|^{2}=\|\hat{x}_{\alpha_{k}}-x^{\dagger}\|^{2}+\|\hat{x}_{\alpha_{k}}^{\delta_{k}}-\hat{x}_{\alpha_{k}}\|^{2}+2\left\langle\hat{x}_{\alpha_{k}}-x^{\dagger},\hat{x}_{\alpha_{k}}^{\delta_{k}}-\hat{x}_{\alpha_{k}}\right\rangle,$

we obtain

$\displaystyle\|\hat{x}_{\alpha_{k}}^{\delta_{k}}-x^{\dagger}\|^{2}\geq\|\hat{x}_{\alpha_{k}}^{\delta_{k}}-\hat{x}_{\alpha_{k}}\|^{2}\geq\delta_{k}^{2}\frac{\lambda_{k}/2}{(\alpha_{k}+3\lambda_{k}/2)^{2}}=\frac{\delta_{k}^{2}/\lambda_{k}}{2(\alpha_{k}/\lambda_{k}+3/2)^{2}}.$
(11)

Since $\delta_{k}=\lambda_{k}$, it then follows from (10) and (11) that

$\frac{\delta_{k}}{2(\alpha_{k}/\delta_{k}+3/2)^{2}}\leq\|\hat{x}_{\alpha_{k}}^{\delta_{k}}-x^{\dagger}\|^{2}=o(\delta_{k})$

which implies that $\alpha_{k}/\delta_{k}\to\infty$ as $k\to\infty$.

Finally we show $x^{\dagger}=x^{*}$. Let $C_{1}:=\sup_{x\in
B_{\rho}(x^{\dagger})}\|F^{\prime}(x)\|$ which is finite by Assumption 1. By
the definition of $\alpha_{k}$ and the Euler equation for
$x_{\alpha_{k}}^{\delta_{k}}$ we have

$\displaystyle\|x_{\alpha_{k}}^{\delta_{k}}-x^{*}\|=\frac{1}{\alpha_{k}}\|F^{\prime}(x_{\alpha_{k}}^{\delta_{k}})^{*}(F(x_{\alpha_{k}}^{\delta_{k}})-y^{\delta_{k}})\|\leq\frac{C_{1}}{\alpha_{k}}\|F(x_{\alpha_{k}}^{\delta_{k}})-y^{\delta_{k}}\|\leq\frac{C_{1}\tau\delta_{k}}{\alpha_{k}}\to
0$

as $k\to\infty$. Since (8) also implies
$\|x_{\alpha_{k}}^{\delta_{k}}-x^{\dagger}\|\to 0$ as $k\to\infty$, we
therefore have $x^{\dagger}=x^{*}$. ∎

The saturation result given in Theorem 2 covers not only Rule 1 but also other
possible parameter choice rules including the sequential discrepancy principle
([6]) below.

######  Rule 2 (Sequential discrepancy principle).

Let $\tau>1$, $\alpha_{0}>0$ and $0<\gamma<1$ be given numbers, and let

$\Delta_{\gamma}:=\\{\alpha_{0}\gamma^{j}:j=0,1,\cdots\\}.$

Define $\alpha(\delta,y^{\delta})$ to be the largest number in
$\Delta_{\gamma}$ such that
$\|F(x_{\alpha(\delta,y^{\delta})}^{\delta})-y^{\delta}\|\leq\tau\delta$.

We remark that for nonlinear ill-posed problems there are some issues
concerning Rule 1 which require serious consideration. First, the existence of
$\alpha(\delta,y^{\delta})$ in Rule 1 is not always guaranteed in general.
Furthermore, even if such an $\alpha(\delta,y^{\delta})$ exists, additional
work needs to be executed to find it which can be time-consuming. The
sequential discrepancy principle can resolve these issues automatically.
Indeed, by the definition of $x_{\alpha}^{\delta}$ it is easy to see that
$\lim_{\alpha\to 0}\|F(x_{\alpha}^{\delta})-y^{\delta}\|\leq\delta$ which
implies the existence of $\alpha(\delta,y^{\delta})$ by Rule 2. According to
Theorem 2, under merely Assumption 1, Tikhonov regularization (3) with Rule 2
exhibits the saturation.

Next we will provide a saturation result under the following Lipschitz
continuity of the Fréchet derivative of the forward operator $F$.

######  Assumption 3.

There is $\rho>0$ such that $B_{\rho}(x^{\dagger})\subset\emph{dom}(F)$ and
$F$ is Fréchet differentiable on $B_{\rho}(x^{\dagger})$. Moreover, there is
$L\geq 0$ such that

$\|F^{\prime}(x)-F^{\prime}(z)\|\leq L\|x-z\|,\qquad\forall x,z\in
B_{\rho}(x^{\dagger}).$

We will show that, under Assumption 3, Tikhonov regularization (3) with the
regularization parameter chosen by Rule 1 or Rule 2 does not possess the
convergence rate $o(\delta^{1/2})$ in the worst case scenario even if
$x^{\dagger}-x^{*}\in\mbox{Ran}((A^{*}A)^{\nu})$ with $\nu\geq 1/2$. We need
the following result which is essentially contained in the proof of [3,
Theorem 2.4], see also [2, 12].

######  Lemma 3.

Let Assumption 3 hold and let $x_{\alpha}^{\delta}$ be a solution of (3) with
$x_{\alpha}^{\delta}\in B_{\rho}(x^{\dagger})$. If
$x^{*}-x^{\dagger}=F^{\prime}(x^{\dagger})^{*}u$ for some $u\in Y$ with
$L\|u\|<1$, then

$\|x_{\alpha}^{\delta}-x^{\dagger}\|\leq\frac{\delta+\alpha\|u\|}{\sqrt{\alpha(1-L\|u\|)}}\quad\mbox{and}\quad\|F(x_{\alpha}^{\delta})-y^{\delta}\|\leq\delta+2\alpha\|u\|.$
(12)

######  Theorem 4.

Let Assumption 3 hold and assume, for the operator
$A:=F^{\prime}(x^{\dagger})$, there exists a sequence of positive numbers
$\\{\lambda_{k}\\}\subset\sigma(AA^{*})$ such that
$\lim_{k\to\infty}\lambda_{k}=0$. Consider the Tikhonov regularization (3) and
let $\alpha(\delta,y^{\delta})$ be determined by either Rule 1 or Rule 2 with
$\tau>1$. If $x^{\dagger}-x^{*}=F^{\prime}(x^{\dagger})^{*}u$ for some
$u\in\mathcal{N}(F^{\prime}(x^{\dagger})^{*})^{\perp}$ with $L\|u\|<1$ and if

$\displaystyle\sup\left\\{\|x_{\alpha(\delta,y^{\delta})}^{\delta}-x^{\dagger}\|:\|y^{\delta}-y\|\leq\delta\right\\}=o(\delta^{1/2})$
(13)

as $\delta\to 0$, then $x^{\dagger}=x^{*}$.

###### Proof.

Let $\alpha:=\alpha(\delta,y^{\delta})$. From (13) it follows that
$x_{\alpha}^{\delta}\in B_{\rho}(x^{\dagger})$ for small $\delta>0$. Thus, by
virtue of (2) and Assumption 3 we have

$\displaystyle\|x_{\alpha}^{\delta}-\hat{x}_{\alpha}^{\delta}\|$
$\displaystyle\leq\frac{1}{2\sqrt{\alpha}}\|F(x_{\alpha}^{\delta})-y-A(x_{\alpha}^{\delta}-x^{\dagger})\|+\frac{1}{\alpha}\|F^{\prime}(x_{\alpha}^{\delta})-A\|\|y^{\delta}-F(x_{\alpha}^{\delta})\|$
$\displaystyle\leq\frac{1}{4\sqrt{\alpha}}L\|x_{\alpha}^{\delta}-x^{\dagger}\|^{2}+\frac{1}{\alpha}L\|x_{\alpha}^{\delta}-x^{\dagger}\|\|y^{\delta}-F(x_{\alpha}^{\delta})\|.$

With the help of Lemma 3 we then have

$\displaystyle\|x_{\alpha}^{\delta}-\hat{x}_{\alpha}^{\delta}\|\leq\frac{L(\delta+\alpha\|u\|)}{4\alpha\sqrt{1-L\|u\|}}\|x_{\alpha}^{\delta}-x^{\dagger}\|+\frac{1}{\alpha}L\left(\delta+2\alpha\|u\|\right)\|x_{\alpha}^{\delta}-x^{\dagger}\|.$
(14)

We next claim that

$\displaystyle\alpha\geq\frac{(\tau-1)\gamma\delta}{2\|u\|}.$ (15)

Indeed, if $\alpha$ is determined by Rule 1, then by using the second
inequality in Lemma 3 we have

$\tau\delta=\|F(x_{\alpha}^{\delta})-y^{\delta}\|\leq\delta+2\alpha\|u\|$

which implies that $\delta\leq 2\alpha\|u\|/(\tau-1)$ and hence (15) holds as
$0<\gamma<1$. If $\alpha$ is determined by Rule 2, then by using again the
second inequality in Lemma 3 we can obtain

$\tau\delta<\|F(x_{\alpha/\gamma}^{\delta})-y^{\delta}\|\leq\delta+\frac{2\alpha}{\gamma}\|u\|$

which implies again (15). Therefore, it follows from (14) and (15) that

$\displaystyle\|x_{\alpha}^{\delta}-\hat{x}_{\alpha}^{\delta}\|\leq
C_{0}L\|u\|\|x_{\alpha}^{\delta}-x^{\dagger}\|,$

where

$C_{0}:=2+\frac{2}{(\tau-1)\gamma}+\frac{2+(\tau-1)\gamma}{4(\tau-1)\gamma\sqrt{1-L\|u\|}}.$

Consequently

$\|\hat{x}_{\alpha}^{\delta}-x^{\dagger}\|\leq(1+C_{0}L\|u\|)\|x_{\alpha}^{\delta}-x^{\dagger}\|.$

and thus we may use (13) to conclude

$\displaystyle\sup\left\\{\|\hat{x}_{\alpha}^{\delta}-x^{\dagger}\|:\|y^{\delta}-y\|\leq\delta\right\\}=o(\delta^{1/2}).$

Now we can follow the same argument in the proof of Theorem 2 to conclude
$x^{\dagger}=x^{*}$. ∎

## References

  * [1] H. W. Engl and H. Gfrerer, A posteriori parameter choice for general regularization methods for solving linear ill-posed problems, Appl. Numer. Math., 4 (1988), pp. 395–417.
  * [2] H. W. Engl, M. Hanke and A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996.
  * [3] H. W. Engl, K. Kunisch and A. Neubauer, Convergence rates of Tikhonov regularization of nonlinear ill-posed problems, Inverse Problems, 5 (1989), 523–540.
  * [4] C. W. Groetsch, Comments on Morozov’s Discrepancy Principle, In: Hämmerlin G., Hoffmann K. H. (eds.), Improperly posed problems and their numerical treatment, pp. 97–104. Basel: Birkhäuser 1983.
  * [5] C. W. Groetsch, The Theory of Tikhonov Regularization for Fredholm Equations of the First Kind, Pitman, Boston, 1984.
  * [6] B. Hofmann and P. Mathé, Parameter choice in Banach space regularization under variational inequalities, Inverse Problems, 28 (2012), no. 10, 104006, 17 pp
  * [7] Q. Jin, Applications of the modified discrepancy principle to Tikhonov regularization of nonlinear ill-posed problems, SIAM J. Numer. Anal., 36 (1999), no. 2, 475–490.
  * [8] Q. Jin and Z. Hou, On an a posteriori parameter choice strategy for Tikhonov regularization of nonlinear ill-posed problems, Numer. Math., 83 (1999), 139–159.
  * [9] A. Neubauer, On converse and saturation results for Tikhonov regularization of linear ill-posed problems, SIAM J. Numer. Anal., 34 (1997), no. 2, 517–527.
  * [10] O. Scherzer, The use of Morozov’s discrepancy principle for Tikhonov regularization for solving nonlinear ill-posed problems, Computing, 51 (1993), no. 1, 45–60.
  * [11] O. Scherzer, H. W. Engl and K. Kunisch, Optimal a posteriori parameter choice for Tikhonov regularization for solving nonlinear ill-posed problems, SIAM J. Numer. Anal., 30 (1993), 1796–1838.
  * [12] U. Tautenhahn and Q. Jin, Tikhonov regularization and a posteriori rules for solving nonlinear ill posed problems, Inverse Problems, 19 (2003), no. 1, 1–21.


# Pronunciation recognition of English phonemes /@/, /æ/, /A:/ and /2/ using
Formants and Mel Frequency Cepstral Coefficients

Keith Y. Patarroyo and Vladimir Vargas-Calderón Physics Department. National
University of Colombia.
Corresponding author e-mail<EMAIL_ADDRESS>
The Vocal Joystick Vowel Corpus, by Washington University, was used to study
monophthongs pronounced by native English speakers. The objective of this
study was to quantitatively measure the extent at which speech recognition
methods can distinguish between similar sounding vowels. In particular, the
phonemes /@/, /æ/, /A:/ and /2/ were analysed. 748 sound files from the corpus
were used and subjected to Linear Predictive Coding (LPC) to compute their
formants, and to Mel Frequency Cepstral Coefficients (MFCC) algorithm, to
compute the cepstral coefficients. A Decision Tree Classifier was used to
build a predictive model that learnt the patterns of the two first formants
measured in the data set, as well as the patterns of the 13 cepstral
coefficients. An accuracy of 70% was achieved using formants for the mentioned
phonemes. For the MFCC analysis an accuracy of 52 % was achieved and an
accuracy of 71% when /@/ was ignored. The results obtained show that the
studied algorithms are far from mimicking the ability of distinguishing subtle
differences in sounds like human hearing does.

Keywords: sound processing, formants, Mel frequency cepstral coefficients,
pronunciation recognition, machine learning.

##  1 Introduction

Throughout computing history, scientists have developed a vast amount of
theories and algorithms for speech recognition that are widely known (e.g. see
the work by Lee (1988)). Many of them are motivated by using some of the
principles of the human ear operation (Davis & Mermelstein, 1980). In recent
years, the blooming of high-speed processing computers has allowed us to start
using deep machine learning to improve the efficacy of our speech recognizers
(Baker et al., 2009). The techniques that are currently used for speech
recognizers yield high prediction rates (Hinton et al., 2012). The most common
and successful technique is the implementation of multi-layered neural
networks (Zegers, 1998; Wellekens, 1998). The performance of such techniques
seeded the question of implementing recognizers as objective evaluators of
speech ability in humans.

This takes great relevance in the field of pronunciation teaching, not only
because it provides the teachers (especially teachers who are non-native
speakers of the taught language) a tool for assessing objectively the
pronunciation of their students (Neri, Mich, Gerosa, & Giuliani, 2008), and of
themselves, but also because if a student has access to such a tool, he or she
could learn pronunciation autonomously (Hinks, 2003). Unquestionably, the main
goal in pronunciation teaching is to improve the ability of a language learner
to use their vocal tract to produce sounds that are recognized as native
sounds by native speakers of that language. Although this remains a challenge,
particularly in the early stages of learning a language (Mirzaei, Gowhary,
Azizifar, & Esmaeili, 2015), it also provides a means by which students can
improve the learning and retention of grammatical structures (Martin &
Jackson, 2016). Also, since different languages have different sound systems,
it is normal for language learners to struggle learning them.

For instance, the phonemes /@/, /æ/, /A:/ and /2/ found in the English
language are troublesome for many English learners because of their
similarity. For example, Spanish and Italian speakers tend to mix these
phonemes into a single one: /a/. There are cases of learners from other
languages, such as Azerbaijani, in which this confusion has been studied
(Ghaffarvand Mokari & Werner, 2016).

This article aims to study the ability of two of the most historically
important speech recognition tools to differentiate between these phonemes,
i.e. how good tools would they be in evaluating the correct pronunciation of
these phonemes. These tools are the formants (computed with Linear Predictive
Coding), attractive for its simplicity in vowel pronunciation recognition, and
the Mel Frequency Cepstral Coefficients, attractive for its computational
speed in continuous speech recognition. Both tools were tested with the Vocal
Joystick Vowel Corpus, from the Washington University.

##  2 Theoretical Framework

###  2.1 Sound Formation

The frequency range of the pressure waves in the air that compose the audible
sounds by humans goes from 20$\mathrm{Hz}$ to 20$\mathrm{kHz}$ (Rosen &
Howell, 2011). The process by which sound is produced by humans is the
following (for a complete understanding of sound production by humans see
reference (O'Grady & Dobrovolosky, 1997)).

The diaphragm and intercostal contractions make the lungs generate a flow of
air from the chest to the mouth. This air first passes through the larynx,
wherein the vocal cords are. These are muscles that can be geometrically
distributed in several ways, each of which is excited by the passing air
creating modes of vibration or glottal states that result in sound. The
pharynx, the oral cavity and the nasal cavity are filters and resonators of
the aforementioned sounds. Also, the tongue and lips allow us to rapidly
articulate and change the shape of the vocal cavity in order to filter some
frequencies.

In this study, we are interested in vowels, which are voiced glottal states
produced with little obstruction of the vocal tract. This means having the
lips open and also the tongue without contact with the palate.

###  2.2 Human Hearing Detection Principles

It is also convenient to mention some of the principles that the human ear
uses to detect sound signals. The human ear analyses pressure variations in
the air into different frequencies. Roughly, the human ear does a Fourier
transform of the pressure signal $P(t)$, where $t$ is the time, and transmits
$|P(\omega)|^{2}$, where $\omega$ is the frequency, to the brain (Sethna,
2006).

However, this initial model falls short for many analyses. Some of the reasons
are (Lyons, n.d.): the tonal information is changing as a word or tune
progresses; the difference between two closely spaced frequencies is hard to
recognize for humans, especially at high frequencies; and humans hear loudness
on a non-linear scale. All of these factors propose a significant challenge to
mimic the human hearing system.

Davis and Mermelstein (1980) developed a method (Mel Frequency Cepstral
Coefficients) to mimic this process, and has been the state of the art in the
field of speech recognition since then.

###  2.3 Formants and Linear Predictive Coding (LPC)

The vocal tract can be thought of as a vibrating cavity, or resonator, whose
geometry (morphology) is constantly changing in continuous speech. This change
in the geometry allows the superposition of different modes of vibration, and
thereby different sounds.

Experiments have shown that for English vowels there are some characteristic
frequencies of vibration called formants (Hillenbrand, Getty, Wheeler, &
Clark, 1994; Hunter & Kebede, 2012; Deterding, 2006), that correspond to
maxima of vibration, i.e. where the acoustic energy is focused. The first
formant is roughly located from 0 to 1$\mathrm{kHz}$, while the second formant
is located from 1$\mathrm{kHz}$ to 2$\mathrm{kHz}$, and so on (these frequency
boundaries are not rigid, because there might be some second formants below
1$\mathrm{kHz}$, as seen in the references). In the studies mentioned only the
two first formants are taken into account for characterising each vowel.
However, there are studies like (Prica & Ilić, 2010) in which 3 formants are
used to characterise the Serbian vowels in order to improve accuracy when
performing classification of vowel sounds.

Formants can be found by performing an acoustic power spectrum of a sound
signal (or a spectrograph of the signal), and identifying the peaks in the
spectrograph. This can be done by computing the envelope of the signal's
frequency spectrum. Since the envelope contains information about the energy
peaks, the formants can be determined. To predict the envelope of a signal
$s[n(\Delta t)]$ sampled in discrete time steps $n=0,1,2,\ldots$, the Linear
Predictive Coding method was introduced (Deng & O'Shaughnessy, 2003). The goal
of LPC is to predict $s[n]$ with a linear combination of
$s[n-1],s[n-2],\ldots,s[n-M]$, where $M$ is the integer that sets the number
of predicting signal samplings.

###  2.4 Mel Frequency Cepstral Coefficients (MFCC)

This method aims to mimic the procedure performed by the human hearing system
to decode a sound signal. It can be summarised in the steps shown in table 1
(Lyons, n.d.).

Table 1: Steps of the MFCC method compared to the sound analysis of the human body. Step |  Human Detection |  Mel Frequency Coefficients
---|---|---
1. |  The hearing sense is sensible to the tonal information change as a word or tune progresses. |  Frame the signal into short frames to account that on short time scales the audio signal does not change much.
2. |  Different frequencies are detected by vibrations at different spots of the cochlea depending on the frequency of the incoming sounds. |  For each frame calculate the periodogram estimate of the power spectrum.
3. |  The cochlea cannot discern the difference between two closely spaced frequencies. |  Apply the Mel filterbank to the power spectra, sum the energy in each filter.
4. |  Loudness is detected on a non-linear scale, where differences in small frequencies are more considerable. |  Take the logarithm of all filterbank energies.
5. |  ——————————— |  Take the DCT (Discrete Cosine Transform) of the log filterbank energies.
6. |  Identifies the linguistic content and discards all the other stuff (background noise, emotion) |  Keep DCT coefficients 2-13, discard the rest.

###  2.5 Decision Tree Classifiers (DTC)

Classifiers are computer tools that learn patterns from labelled data, so that
when a new sample is presented to the computer, it will classify the sample
into one of the different labels. In other words, consider a data set
$D=\\{(d_{i},l_{i})\\}$, where $d_{i}$ is a vector of $N$ components
(features). The vector $d_{i}$ is also called a sample. $l_{i}$ is a label.
The classifier is a function $f(d):F\to L$, where $F$ is the vector space of
the samples and $L$ the space of the labels.

To see how DTCs work, consider the following example. Suppose that a data set
consists of samples in one dimension: each sample is the salary of every
single person in Colombia. The labels for this data set are the social
stratum. We might have people that earn a lot of money, but they live in
stratum 2, or vice versa. However, these cases are strange and in general
there is a structure that allows predicting the label (social stratum) from
the samples. The DTC will try to learn the salary ranges for each single
social stratum, so that when a new sample is presented to the model, it will
classify it in a specific label. The boundaries of the salary ranges are
random at first, but as the learning process advances, the boundaries become
more accurate and clear.

To test how well the DTC is able to predict labels, a data set $D^{\prime}$
called the testing data is normally built. $D^{\prime}$ consists of $M$
samples with known labels. Only the samples, and not the labels, are presented
to the DTC. The percentage of correctly predicted labels will tell the
accuracy or effectiveness of the DTC. For a more thorough explanation, refer
to (Tan, Steinbach & Kumar, 2005).

##  3 Methodology

The pipeline of this study (see Figure 1) was to extract and prepare audio
files from the Vocal Joystick Vowel Corpus. Then, formants and cepstral
coefficients were computed for each file. The set of files with their
corresponding formants and cepstral coefficients was split into training and
testing data. A DTC was trained with the training data, and tested with the
testing data, yielding an accuracy percentage. The characteristics of each
step are explained in this section.

Figure 1: Diagram of the method used to measure the performance of formants
and MFCC in pronunciation prediction.

###  3.1 Corpus and Audio Processing

In this study we used the Vocal Joystick Vowel Corpus (Bilmes, Wright, Xiao,
Malkin Kilanski, 2006). This project selected a group of nine monophthongs and
12 vowel- to-vowel transitions. Sounds with different duration, amplitude and
intonation were recorded for each vowel:

  * •

Duration: short (1 second), long (2 seconds) and nudge (very short repetitions
of the same vowel).

  * •

Amplitude: quiet, normal, loud, quiet to loud and loud to quiet.

  * •

Intonation: level, rising and falling.

From this corpus, we considered all the sound files available in the corpus
corresponding to the phonemes /@/, /æ/, /A:/ and /2/ whose duration was either
short or long, whose amplitude was quite, normal or loud, and whose intonation
was in a single level. A total of 784 files satisfied these conditions.

###  3.2 Data Preparation

An example of a raw data file from the corpus is shown in Figure 2a. It can be
seen that there are silent parts and the amplitude is not normalised. Also,
even though we selected a single level of intonation, it is clear that the
amplitude of the sound decreases as time progresses. Therefore, for each file
we deleted the silent part. Then, a section of duration 40ms was selected so
that the signal does not change considerably and is approximately periodic
(Figure 2b). Finally, a Hamming window was applied to the resulting sound
signal with the purpose of smoothing the Fourier Transform used in the MFCC
extraction (Figure 2c).

Figure 2: Preparation example of a sound signal corresponding to the /2/
phoneme. a) Raw signal. b) Silence deletion, amplitude normalisation and
selection of a 40ms interval. c) Application of a Hamming window.

###  3.3 Evaluation of Formants and MFCC Performance with a DTC

Then, the LPC Python implementation by Danilo Bellini was used to compute the
formants (Bellini, 2016). To do the calculation of the MFCC, the Python
implementation by James Lyons (Lyons, 2016) was used. Finally, the set of
computed formants and the set of computed cepstral coefficients were
partitioned into two groups each. One of the groups contained two thirds of
the data, and the other one contained a third. The former was the training
data, and the later was the testing data. To exemplify this procedure, denote
$S=\\{(s,v)\\}$ as the set of the sound signal files, each of which is
labelled by a vowel $v$. Let $F(s)$ be a function that computes the formants
$F$ of a sound signal $s$. Let $X=\\{(F(s),v)\,\forall s\in S\\}$ be the data
set of formants, labelled by their respective vowel. Let $X_{tr}\subset X$ be
the set of training data and $X_{te}\subset X$ be the set of testing data. The
DTC learns the ranges corresponding to each vowel in the formants space by
creating a function $DTC(F(s))$ that takes the values of the labels, i.e. the
phonemes /@/, /æ/, /A:/ and /2/. For instance, if after the learning process
is finished, a sample $F(s)$, labelled by the phoneme /@/, is presented to the
DTC, and $DTC(F(s))=$ /æ/, then the DTC failed to predict the correct phoneme.
If instead $DTC(F(s))=$ /@/ then the DTC was successful in the prediction. The
accuracy is thus defined as the percentage of correctly predicted phonemes
from the testing data. We used the accuracy to test the quality of prediction
by the DTC using formants and Mel coefficients separately. This gives a
measure of how well could formants and cepstral coefficients help in the
objective evaluation of pronunciation.

##  4 Results, Analysis and Discussion

###  4.1 Formants

From the 784 sound files, 142 were discarded for this analysis and for the one
corresponding to the MFCC. The reason to discard these files was that some
recorded sounds contained parts in which the signal was not periodic due to a
trembling voice from the speaker. This could cause the LPC and MFCC algorithms
to be affected. The criterion used to discard the data is now explained.
Denote the recordings or sound signals that belong to vowel $v$ as $(s,v)$.
Also, denote the mean value of the $i$-th formant ($i=1,2$) and the standard
deviation of the $i$-th formant, for a vowel $v$, as
$(\mu_{i},v),(\sigma_{i},v)$, respectively. The remaining 642 sound signals
satisfied the condition

$\displaystyle|F_{i}(s)-\mu_{i}|<1.5\sigma_{i},$ (1)

for each vowel $v$, where $F_{i}$ refers to the $i$-th formant. The computed
formants for these sound signals are shown in figure 3.

Figure 3: Frequencies in Hertz of the first two formants for the studied
phonemes. aa is /A:/, ax is /@/, ae is /æ/ and ah is /2/.

The DTC trained with the formants yielded a 70% accuracy. Although this
success rate of the predictive model is a very low rate, it can be considered
as good taking into account that only two features were used to predict the
data: the first two formants. This means that with little information about
four very similar sounding phonemes, we were able to predict, with an accuracy
of 70%, the vowel of a sound taken at random from the recordings obtained from
the corpus.

Furthermore, Figure 3 shows a non-Gaussian distribution in either formant for
every phoneme. This is supported by observations in which it was seen that
formants do not only depend on the geometry of the vocal tract of each person,
but also are statistically dependent on the age, as shown in (Hawkings &
Midgley, 2005); as well as on gender, as shown in (Hillenbrand, Getty, Clark,
& Wheeler, 1995). The formants measured are shown in table 2 and compared with
measurements from other studies.

Table 2: Formants measurement in Hertz. Our measurements are compared with Hawkings et. al., Hunter et. al., Deterding. Study | Our | Hawkings et. al. | Hunter et. al. | Deterding
---|---|---|---|---
Phoneme | $F_{1}$ | $F_{2}$ | $F_{1}$ | $F_{2}$ | $F_{1}$ | $F_{2}$ | $F_{1}$ | $F_{2}$
/@/ | 631 | 1049 | 496 | 833 | 643 | 1019 | 625 | 973
/æ/ | 720 | 1644 | 696 | 1574 | 667 | 1565 | 748 | 1360
/A:/ | 573 | 1311 | 608 | 1062 | 680 | 1193 | 757 | 1211
/2/ | 693 | 1182 | 629 | 1160 | 661 | 1296 | 724 | 1282

Studies shown in Table 2 show that formants vary a lot with age, gender and
geography. Therefore, we see that identifying a phoneme by its formants is not
trivial, but as we showed, can be done with an estimated accuracy of 70%.
Despite the low accuracy, formants keep being used for vowel pronunciation
research, as in (Rasilo & Räsänen, 2017), in which the process of infants
learning their native language was simulated through a Leaning Virtual Infant
that received sounds from its caregivers. At first the Learning Virtual Infant
babbled randomly, but as the learning process advanced, it began to compute
formants to identify similar sounding vowels and to improve its pronunciation.

###  4.2 MFCC

13 MFCC were obtained for every sound signal. In order to visualise the MFCC
of every sound signal, a dimension reduction from 13 to 2 (i.e. from the MFCC
space to a plane) was performed with Principal Component Analysis (PCA). This
method projects the data from the high-dimensional space onto a plane that
preserves the maximum possible variance of the data. The distribution of the
MFCC for the different phonemes on the plane computed with PCA is shown in
Figure 4.

Figure 4: PCA of the MFCC. a) shows the four studied phonemes and b) shows
only /@/, /æ/ and /2/.

As it can be seen from Figure 4a, the MFCC of the phonemes seem mixed. In
particular, note that the points corresponding to the phoneme /A:/ are mixed
with the ones corresponding to the phoneme /2/. This can also be seen in the
formants figure. Ignoring the /A:/ points produces Figure 4b, in which the
clusters of data can be identified.

Partitioning the data shown in Figure 4a in training and testing data, an
accuracy of 56% was accomplished with the DTC. While if the data shown in
Figure 4b was used, an accuracy of 71% was measured. These very low accuracies
show that the MFCC method lacks sensibility to subtle differences between the
studied sound signals.

But, why are the MFCC so good at speech recognition then? The first thing to
mention is that the MFCC are normally used to train models that are more
complex than a DTC, like neural networks. The goal in speech recognition is to
understand words, instead of single sounds. To check the accuracy of the MFCC
at recognising different sounding sounds, we proceeded to analyse the phonemes
/2/, /e/, /u/ and /o/. The PCA results are shown in Figure 5.

Figure 5: PCA of the MFCC. a) shows the phonemes /2/, /e/, /u/ and /o/. b)
shows only /2/, /e/ and /u/. ee is /e/, uu is /u/ and oo is /o/.

Using the data shown in Figure 5a, the accuracy was measured in 73%. Clearly
the MFCC corresponding to the phoneme /o/ were mixed with the ones
corresponding to the phoneme /u/. To avoid this noise, if /o/ is removed, an
accuracy of 88% was reached. This shows that MFCC are a reasonable method to
recognise different phonemes.

##  5 Conclusions and Perspectives

The Vocal Joystick Corpus was used to build a data set of sounds corresponding
to the phonemes /@/, /æ/, /A:/ and /2/. From each of the sound signals, both
formants and MFCC were computed. These were used to train a DTC that acted as
a classifier. Using formants, an accuracy of 70% was achieved. Using MFCC, an
acuraccy of 52% was measured. MFCC's performance was much better when the DTC
was trained with the coefficients computed from the phonemes /2/, /e/ and /u/,
reaching 88%. However, these speech recognition tools fail to correctly
predict the phonemes of similar sounding signals.

As a response to this deficiency, computer and linguistics scientists work to
build computational tools that learn subtle differences between sounds in
order to assess a person’s pronunciation correctly. For instance, a research
in the Chinese language showed results that improved significantly the
previous recognition models in detecting mispronunciation of phonemes (Wei,
Hu, Hu, & Wang, 2009). However, the advance of technology to support
pronunciation learning must be accompanied by designed learning programs that
help students to learn autonomously. More people need to be encouraged to work
in this interdisciplinary field that is pioneering and improving language
teaching.

## References

  * [1] Baker, J., Li Deng, Glass, J., Khudanpur, S., Chin-hui Lee, Morgan, N., & O'Shaughnessy, D. (2009). Developments and directions in speech recognition and understanding, Part 1 [DSP Education]. IEEE Signal Processing Magazine, 26(3), 75-80. http://dx.doi.org/10.1109/msp.2009.932166
  * [2] Bellini, D. (n.d.). Expressive Digital Signal Processing (DSP) package for Python, 2016. Retrieved November 16, 2016, from https://github.com/danilobellini/audiolazy
  * [3] Bilmes, J., Wright, R., Xiao, L., Malkin, J. & Kilanski, K. 2006. http://melodi.ee.washington.edu/vj/index.html
  * [4] Davis, S. & Mermelstein, P. (1980). Comparison of parametric representations for monosyllabic word recognition in continuously spoken sentences. IEEE Transactions On Acoustics, Speech, And Signal Processing, 28(4), 357-366. http://dx.doi.org/10.1109/tassp.1980.1163420
  * [5] Deng, L. & O'Shaughnessy, D. (2003). Speech processing (1st ed., pp. 19-26). New York: Marcel Dekker.
  * [6] Deterding, D. (2006). The North Wind versus a Wolf: short texts for the description and measurement of English pronunciation. Journal Of The International Phonetic Association, 36(02), 187. http://dx.doi.org/10.1017/s0025100306002544
  * [7] Ghaffarvand Mokari, P. & Werner, S. (2016). Perceptual assimilation predicts acquisition of foreign language sounds: The case of Azerbaijani learners’ production and perception of Standard Southern British English vowels. Lingua, 185, 81-95. http://dx.doi.org/10.1016/j.lingua.2016.07.008
  * [8] Hawkins, S. & Midgley, J. (2005). Formant frequencies of RP monophthongs in four age groups of speakers. Journal Of The International Phonetic Association, 35(02), 183. http://dx.doi.org/10.1017/s0025100305002124
  * [9] Hillenbrand, J., Getty, L., Clark, M., & Wheeler, K. (1995). Acoustic characteristics of American English vowels. The Journal Of The Acoustical Society Of America, 97(5), 3099. http://dx.doi.org/10.1121/1.411872
  * [10] Hincks, R. (2003). Speech technologies for pronunciation feedback and evaluation. Recall, 15(01). http://dx.doi.org/10.1017/s0958344003000211
  * [11] Hinton, G., Deng, L., Yu, D., Dahl, G., Mohamed, A., & Jailty, N. et al. (2012). Deep Neural Networks for Acoustic Modeling in Speech Recognition. Retrieved from http://static.googleusercontent.com/media/research.google.com/en//pubs/archive/38131.pdf
  * [12] Hunter, G. & Kebede, H. (2012). Formant frequencies of British English vowels produced by native speakers of Farsi. Société Française d'Acoustique. Retrieved from https://hal.archives-ouvertes.fr/hal-00810580/document
  * [13] Lee, K. (1988). On large-vocabulary speaker-independent continuous speech recognition. Speech Communication, 7(4), 375-379. http://dx.doi.org/10.1016/0167-6393(88)90053-2
  * [14] Lyons, J. Mel Frequency Cepstral Coefficient (MFCC) tutorial, Practical Cryptography. Retrieved October 15, 2016, from http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/
  * [15] Martin, I. & Jackson, C. (2016). Pronunciation Training Facilitates the Learning and Retention of L2 Grammatical Structures. Foreign Language Annals, 49(4), 658-676. http://dx.doi.org/10.1111/flan.12224
  * [16] Mirzaei, K., Gowhary, H., Azizifar, A., & Esmaeili, Z. (2015). Comparing the Phonological Performance of Kurdish and Persian EFL Learners in Pronunciation of English Vowels. Procedia - Social And Behavioral Sciences, 199, 387-393. http://dx.doi.org/10.1016/j.sbspro.2015.07.523
  * [17] Neri, A., Mich, O., Gerosa, M., & Giuliani, D. (2008). The effectiveness of computer assisted pronunciation training for foreign language learning by children. Computer Assisted Language Learning, 21(5), 393-408. http://dx.doi.org/10.1080/09588220802447651
  * [18] Oppenheim, J. & Magnasco, M. (2013). Human Time-Frequency Acuity Beats the Fourier Uncertainty Principle. Physical Review Letters, 110(4). http://dx.doi.org/10.1103/physrevlett.110.044301
  * [19] O'Grady, W. & Dobrovolosky, M. (1997). Contemporary Linguistics (3rd ed., pp. 40-87). Boston: Bedford/St. Martin's.
  * [20] Prica, B. & Ilić, S. (2010). Recognition of vowels in continuous speech by using formants. Facta Universitatis - Series: Electronics And Energetics, 23(3), 379-393. http://dx.doi.org/10.2298/fuee1003379p
  * [21] Rasilo, H. & Räsänen, O. (2017). An online model for vowel imitation learning. Speech Communication, 86, 1-23. http://dx.doi.org/10.1016/j.specom.2016.10.010
  * [22] Rosen, S. & Howell, P. (2011). Signals and systems for speech and hearing (2nd ed., p. 163). London: Academic Press.
  * [23] Sethna, J. (2006). Statistical mechanics: Entropy, Order Parameters and Complexity (1st ed., p. 299). Oxford: Oxford University Press.
  * [24] Tan, P., Steinbach, M., & Kumar, V. (2005). Introduction to Data Mining (1st ed., p. 150). Boston: Pearson.
  * [25] Wei, S., Hu, G., Hu, Y., & Wang, R. (2009). A new method for mispronunciation detection using Support Vector Machine based on Pronunciation Space Models. Speech Communication, 51(10), 896-905. http://dx.doi.org/10.1016/j.specom.2009.03.004
  * [26] Wellekens, C. (1998). Introduction to Speech Recognition Using Neural Networks. Bruges: ESANN. Retrieved from https://www.elen.ucl.ac.be/Proceedings/esann/esannpdf/es1998-456.pdf
  * [27] Zegers, P. (1998). Speech Recognition Using Neural Networks (Master of Science with a major in Electrical Engineering). University of Arizona.