Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

55, 4, pp. 1269-1278, Warsaw 2017
DOI: 10.15632/jtam-pl.55.4.1269

COMPARISON OF BAYESIAN AND OTHER APPROACHES TO THE

ESTIMATION OF FATIGUE CRACK GROWTH RATE FROM

2D TEXTURAL FEATURES

Matej Mojzeš, Jaroḿir Kukal

Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Department of Software

Engineering, Praque, Czech Republic

e-mail: mojzemat@fjfi.cvut.cz

Hynek Lauschmann

Czech Technical University in Prague, Faculty of Nuclear Sciences and Physical Engineering, Department of Materials,

Praque, Czech Republic

The fatigue crack growth rate can be explained using features of the surface of a structure.
Among other methods, linear regression can be used to explain crack growth velocity. Non-
linear transformations of fracture surface texture features may be useful as explanatory
variables. Nonetheless, the number of derived explanatory variables increases very quickly,
and it is very important to select only fewof thebest performingones andpreventoverfitting
at the same time. To perform selection of the explanatory variables, it is necessary to assess
quality of the given sub-model.We use fractographic data to study performance of different
information criteria and statistical tests as means of the sub-model quality measurement.
Furthermore, to address overfitting,weprovide recommendationsbasedona cross-validation
analysis. Among other conclusions, we suggest the Bayesian Information Criterion, which
favours sub-models fitting the data considerably well and does not lose the capability to
generalize at the same time.

Keywords: quantitative fractography, optimization, heuristic, linear regression, sub-model
selection

1. Introduction

One of the tasks of quantitative fractography consists of modelling the relation between the
fatigue crack growth rate (CGR) and textural features of images of fatigue fracture surfaces as
explained in (Lauschmann and Siska, 2012; Nadbal et al., 2008; Lauschmann et al., 2006). For
this purpose, e.g. amultilinear regressionmodel (Sekeresova andLauschmann, 2008;Kunz et al.,
2010) or a neural network may be used. Of these two, the neural network allows us to analyze
and describe the structure of the obtained modeland better imagine the textural subset which
is mutually related with the CGR (Lauschmann and Goldsmith, 2009).
The parameters of the respective regressionmodelmay be estimated using the least squares

method. However, in real-world applications, the basic linear model is not flexible enough to
fit the data. This can be solved by adding terms defined by non-linear functions of the basic
features, e.g. logarithm, square root, etc. However, adding such features is soon limited by the
given number of images.
According toLauschmannandGoldsmith (2009), onepossiblewayaround this limitation is a

two-phase stepwise regressionwith the first stage beingbottom-up stepwise regression beginning
with a constant model and terminating at a given overfitting level p0. In each iteration, a new
explanatory variable is included – the one which maximally decreases the sum of squares of
residui. The second stage is top-down stepwise regression beginning with the final sub-model



1270 M.Mojzeš et al.

from the first stage and terminating at the given final overfitting level pF . In this procedure, an
explanatory variable is selected for the elimination via theWald test on a selected critical level.
While keeping in mind the relevant motivation for this problem, we suggest that instead

of the stepwise regression, an alternative statistical approach based on the method of sub-
-model multiple testing may provide better results (Mojzeš et al., 2012). There is a vast set of
possible criteria that evaluate the quality of a given sub-model and are to beminimized. Further
in the paper, we elaborate on the selection and assessment of some of the criteria. They are
are interesting in the fractographic context, but may be applied generally to multi-parametric
recognition as well.

2. Material and methods

2.1. Linear model

Denote by vj the crack growth rate assigned to the j-th image of the fracture surface, and
by fuj the set of image features. Themultilinearmodel in its basic form is based on the formula

log10vj ≈ c0+
∑

u

cufuj (2.1)

Parameters cu can be estimated by the least squares method. Since the linear model is not
flexible enough to fit the data, we may add different non-linear functions of basic features and,
therefore, modify the model to the following form

log10vj ≈ c0+
∑

q

cqhq (2.2)

where h are selected from an extended set of features containing the features fu and a selection
of their basic non-linear functions, e.g.

{hv}⊂P∪Q∪R (2.3)
where

P= {fu} Q= {log10fu,f−1u ,f1/2u ,f2u}, R= {Fuv,F−1uv ,F1/2uv ,F2uv |u>v}
(2.4)

where Fuv = fufv.
The next task will consist of defining a specific methodology on how to select and assess a

distinct combination of explanatory variables, i.e. how to select the best sub-model from the
extended feature set.

2.2. Sub-model selection

A sub-model should be considered a nested subset of the full model consisting of all the
variables from the extended feature set. There are two extreme cases: the first is the full model
and the second one corresponds to the constant model.
Letn∈ N be the length of v (i.e. thenumberof observations),m∈ N the extended feature set

cardinality and k∈{0,1, . . . ,m} the number of explanatory extended feature set variables used
in the sub-model. Also, let c= [c0,c1, . . . ,ck] be the vector representing sub-model coefficients
calculated solvingEq. (2.2).Thisvector canbedivided intocred = [c1,c2, . . . ,ck] as its significant
part and c0 = [c0] as the constant term coefficient.
Furthermore, wemay express the sum of squares for the optimum c of a given sub-model as

SSQ and SSQ0 as the sum of squares for c0. Lastly, we will use the sub-model error defined as

s
2
e =

SSQ

n−k−1
(2.5)



Comparison of Bayesian and other approaches to the estimation of fatigue crack... 1271

2.3. Data description

The methods developed in this paper will be applied to fatigue fracture surfaces of three
laboratory specimens of the heat-resistant steel P92. Compact tension specimens (Fig. 1) (Lau-
schmann et al., 2011) were loaded in air at 20◦C by constant sinusoidal cycles with parameters
of the external force according to Table 1. The loading frequency was lowered in steps from
13Hz to 4Hz in the final stage.

Fig. 1. Compact tension specimen

Table 1.External force parameters

Specimen Fmin [N] Fmax [N]

1 140 3300

2 2000 4800

3 3300 5500

The fatigue crack surfaces were documented using SEMwith magnification 200×, real field
of view was 0.6mm×0.45mm (examples in Fig. 2). The sequence of images was located in the
middleof the crack surface along the sameaxis according towhich the crack lengthwasmeasured
(Fig. 3). The recorded areas were mutually shifted by 0.4mm. The direction of crack growth in
the imageswas bottom-up.Digital representation in 1200×1600 pixels and 256 brightness values
was used. The estimates of CGR were computed from frequently repeated records of the crack
length. The course of the CGR related to the crack length was estimated, and every image was
assigned a value of the CGR pertinent to its middle.

Fig. 2. Examples of SEMpictures of fracture surface; (a) low crack growth rate, (b) high crack
growth rate



1272 M.Mojzeš et al.

Fig. 3. Layout of snaps in crack surface (schematic plot)

For image textural features, energies of a 2D discrete wavelet transformwere taken (Lausch-
mannandGoldsmith, 2009). Decomposition using theType3Daubechieswavelet at 8 levels was
computed by Matlab function wavedec2. The energy is the mean square of wavelet coefficients
for a given level and direction.
The basic sequence of features, x1,x2, . . . ,x24, may be regarded as a set of

H1,V1,D1, . . . ,H8,V8,D8 where Hj,Vj,Dj are wavelet decomposition energies at the j-th le-
vel in the horizontal, vertical and diagonal directions. The vector y represents decimal logarithm
of the crack growth rate y= log10v.
The analysed data consisted of n = 162 observations and a total of 1224 features in the

expanded feature set. It comprises:

• basic linear featuresP, card(P) = 24,
• non-linear transformations of the basic features Q, card(Q)= 96,
• dot product transformations of the basic featuresR, card(R)= 1104,

as stated in (2.3).
To minimize potential numerical errors when working with the data, input data standardi-

zation was implemented as follows

xk =
hk−Eh√
Dh

(2.6)

using Eh and Dh as the mean value and dispersion of the explanatory data.
Last, but not least – apart from the significance of the data, we canmake use also of physical

distribution of the data in the given data set, which is divided randomly into three separate
groups. This will be especially useful when dealing with the cross-validation.

2.4. Selection heuristic

Searching for the best available sub-model is a binary optimization task that can be defined
as minimization of the objective function f:D→ R where

D= {x∈{0,1}m | 0¬x¬1} (2.7)

is thebinarydomain.Here, thebinaryvectorx is directly representingutilization of the extended
feature set, i.e. its components that are equal to “one” are included in the corresponding sub-
-model. Therefore, 0 refers to the constant model, and 1 to the full model.
Furthermore, suppose that we have an acceptable value of the objective function f∗. Then,

we can define a set of solutions, the goal set, as

G= {x∈D | f(x)¬ f∗} (2.8)

where

f
∗ ­min{f(x) | x∈D} (2.9)



Comparison of Bayesian and other approaches to the estimation of fatigue crack... 1273

For that purpose, we may utilize some of the well-known heuristic algorithms. We have chosen
physically motivated Fast Simulated Annealing (FSA) (Kvasnička et al., 2000) with reputable
efficiency in the case of integer optimization tasks. FSA performs mutation on the ring neigh-
bourhood

N(x)= {y∈D |‖y−x‖1 =1} (2.10)

Beginning with k = 0, Tk > 0 and the initial solution vector generated by uniform distribu-
tion x0 ∼ U(D), we perform FSA mutation as a uniformly generated random binary vector
yk ∼U(N(xk)). Using ηk ∼U([−1,+1]), we set

xk+1 =











yk for f(yk)<f(xk)+Tk tan
πη

2

xk for f(yk)­ f(xk)+Tk tan
πη

2

(2.11)

until a solution from the goal set is found, or the pre-defined number of objective function
evaluations is exhausted. The cooling strategy is represented by a non-increasing sequence of
positive temperatures Tk.

We have been inspired by the increased efficiency of hybrid heuristics in the case of combina-
tion of a differential evolution and the steepest descent (Tvrd́ik J and Křivý, 2011). And since
the previously defined set of optimization problems has many local minima, we have enhanced
the FSA algorithm by a hybrid part – the steepest descent, whichmay increase the probability
of reaching the global optimum.

In our approach to hybrid heuristic optimization, instead of f(x) optimization, we optimize
g(x) = f(h) where x = x0, h = xH are the first and last members of any series {xk}Hk=0
satisfying xj ∈ N(xj−1), f(xj) < f(xj−1) for j = 1, . . . ,H. This means h is the best solution
in terms of the steepest descent heuristic. Before any problem solution vector is evaluated, its
nearest local neighbourhood is iteratively searched for a better solution until no further advance
in terms of the objective function can bemade (or until a pre-definedmaximumnumber of local
evaluations is exceeded).

This way we were able to set a higher temperature T0 and use a more benevolent cooling
strategy. In other words, the algorithm is able to prevent getting stuck in a local minimum and
still not lose the ability to fine-tune a given solution. Thus the FSAperformance, on this specific
task, has been improved.

2.5. Cross-validation

Having the data divided into three distinct groups allows us to perform strong cross-
-validation to assess how the results of a specific criterion will generalize to an independent
data set.

We will perform the optimization on two out of three groups (training groups) and validate
theanalysis on the remainingthirdgroup(verification group).To improve theoverall consistency,
multiple rounds of cross-validation will be performed using different permutations of the data
sets, and the verification results will be aggregated over the rounds.

As the goodness of fit measure we propose to use R as the correlation coefficient between
the original data and the data proposed by the respective sub-model. However, during the
optimization itself, we will use the original objective function based on the minimization of the
CRIT value.



1274 M.Mojzeš et al.

3. Theory of sub-model selection

At this point, we should choose some of many possibilities for testing sub-model quality. We
have selected a few of them in accordance with (Mojzeš et al., 2016) that can be divided in two
sets, based on the concepts they are based on.Thefirst one comprises traditional statistical tests
and the criterion that will reflect the quality of a sub-model will be the logarithm of the pvalue.
On the other hand, the second set contains different statistical information criteria regarding
selection of the model. In the latter case, we are simply minimizing the value of the selected
information criterion.

3.1. Sub-model testing

Here, we test significance of a sub-model described by its cred. Respective hypotheses are
defined as

H0 : cred =0

H1 : cred 6=0

and theR-square as well as theWald test are used for their testing.

R-square test

In order to use R2 for the analysis of the sub-model and constant model, according to the
variance analysis (Wooldridge, 2002), we must define the stochastic variable F

F =
SSQ0−SSQ
SSQ

n−k−1
k

(3.1)

Here F belongs to Fk,n−k−1 and the respective pvalue =1−Fk,n−k−1(F).

Wald test

Alternatively, to test sub-model hypotheses via theWald test (Anděl, 1978), the variableZ
has to be used

Z =
1

ks2e
c
T
W
−1
c (3.2)

The matrix W is a result of (XTX)−1 without the first row and column. Finally, Z has the
distribution Fk,n−k−1 and pvalue =1−Fk,n−k−1(Z).
Lastly, for both tests, the resulting sub-model quality criterion can be defined as

CRIT = log10pvalue (3.3)

to beminimized.Due to the fact that values of pvalue mayget very close to zero, it is necessary to
avoid potential numerical issues andexpresspvalue in termsof an incomplete gammadistribution.

3.2. Information criteria

Another approach to the assessment of the sub-model quality is based on statistical infor-
mation criteria. The selected criteria are presented in order from the least stringent one.



Comparison of Bayesian and other approaches to the estimation of fatigue crack... 1275

Wilks Information Criterion

Ralston and Rabinowitz (2001) according to Wilks (1962) recommend searching for a sub-
-model having the minimal error s2e. The corresponding logarithmic form, consistent with the
following criteria, can be defined as

WIC =n lns2e (3.4)

In this basic criterion, k (i.e. the number of explanatory variables included in the sub-model) is
already indirectly penalizing the information quality.

Akaike Information Criterion

Furthermore, an additional penalty for addingmore explanatory variables is included in the
Akaike criterion measuring relative goodness of the sub-model (Akaike, 1974) denoted as

AIC =WIC +2k (3.5)

Bayesian Information Criterion

Bayesian criterion (Schwarz, 1978) prevents overfitting even more by generating a stronger
penalty for extra explanatory variables. Following the existing terminology, the criterion may
be defined as

BIC =WIC +k lnn (3.6)

for n­ 8.
As opposed to the logarithm of pvalue, the final CRIT to be minimized is directly equal to

values of the information criteria.

4. Results and discussion

4.1. Heuristic optimization

In order to compare the results achievedwith the hybridheuristic to the stepwise alternative,
we have implemented a traditional stepwise approach (Mojzeš et al., 2012). Despite using only
125 variables (basic linear and non-linear transformations, P and Q) the heuristic sub-model
optimization approach was superior to the stepwise approach based on the best correlation
coefficient found.

Table 2 aggregates the best results achieved with the current, much more computationally
demanding, data. In multiple runs of the heuristic optimization, we were able to obtain even
better values of the correlation coefficient. Still, as the ultimate target should be the ability
to generalize the independent data set, we will draw final conclusions based on the following
cross-validation.

4.2. Cross-validation

The full data set has been divided into three groups of data, each having 59 (Group I),
53 (Group II) and 50 (Group III) observations. For each permutation of training and verification
groups the hybrid heuristic optimize the sub-model to make the model fit the training data as
well as possible according to the respectivemethod.The same settings and conditions were used
as in the case of the full data set without cross-validation. The detailed results are organized



1276 M.Mojzeš et al.

Table 2.Optimal sub-model quality and features using hybrid heuristic

Method CRIT R kopt
Term

fu f
1/2
u f

2
u f

−1
u log10fu Fuv F

1/2
uv F

2
uv F

−1
uv

F-test −117.92 0.9909 27 0 0 0 0 0 7 4 11 5
Wald test −100.29 0.9839 15 0 0 0 1 0 1 3 7 3
WIC −1099.00 0.9993 88 1 1 4 1 1 14 19 22 25
AIC −927.03 0.9992 82 0 1 2 1 0 20 19 18 21
BIC −626.86 0.9871 20 0 0 1 0 1 1 2 13 2

Table 3.Cross-validation results

Method Training groups CRIT Rtrain Rverify kopt

R2 test I+II −86.72 0.9970 0.5937 34
R2 test II+III −76.70 0.9900 0.5785 13
R2 test I+III −79.29 0.9924 0.9594 22
Wald test I+II −74.29 0.9867 0.8624 10
Wald test II+III −71.34 0.9848 0.6661 5
Wald test I+III −66.24 0.9854 0.9476 13
WIC I+II −736.29 0.9991 0.4471 55
WIC II+III −804.83 0.9997 0.5214 61
WIC I+III −716.25 0.9993 0.8286 63
AIC I+II −619.41 0.9988 0.8292 48
AIC II+III −674.31 0.9996 0.1351 56
AIC I+III −584.22 0.9987 0.9229 48
BIC I+II −479.12 0.9984 0.6896 45
BIC II+III −500.70 0.9995 0.3800 53
BIC I+III −430.91 0.9926 0.9574 21

in Table 3. The most important results are in the column of the correlation coefficient Rverify
which measures the quality of the fit on the verification data set.

These results are aggregated using themean of respectivemethods and, furthermore, expan-
ded by comparing the data composed from distinct verification data sets to the original one in
Table 4. Also, the results of composed cross-validation are depicted in Fig. 4 as log-log plots of
measuredandpredictedCGRs.Thepredicteddata are distinguishedbygroup symbols (×,+,◦).
The results ofWIC and AIC are the worst ones, since the outliers had to be omitted to enable
the plotting. As can be seen, theWald approach and BIC criterion offer the best results.

Table 4.Cross-validation summary

Method Mean R ComposedR

R2 test 0.7105 0.7043

Wald test 0.8254 0.7702

WIC 0.5990 0.5320

AIC 0.6291 0.3281

BIC 0.6757 0.6096



Comparison of Bayesian and other approaches to the estimation of fatigue crack... 1277

Fig. 4. Prediction of CrackGrowth Rate (CGR) using various techniques of sub-model selection
(× representing data in group I, + in group II, and ◦ in group III)

5. Conclusions

The benefits of the solution described above are considerable. Nearly an unlimited set of expla-
natory variables may be offered without any respect to the original number of observations in
a given case. Very good models have been obtained also in previously unsolvable cases with a
very small number of observations.

Of course, the final result is mostly dependent on the approach of sub-model selection. As
it is apparent from the results of cross-validation, and also based on our experience, we are
recommending BIC, Wald test and potentially also R2 test and WIC. Nevertheless, there are
significant differences between these four and, more specifically, we are suggesting:



1278 M.Mojzeš et al.

• BIC as a universal criterion,
• Wald test as a well balanced criterion, similar to BIC,
• R2 test as a legitimate criterion with respect to the variance analysis approach,
• WIC as a criterion that leads to considerable adherence to the data, however, as opposed
to other criteria, lacks the ability of generalization.

Acknowledgement

This paper was created under the support of grant OHK4-165/11CTU in Prague.

References

1. Akaike H., 1974, A new look at the statistical model identification, IEEE Transactions on Auto-
matic Control,AC-19, 716-723

2. Anděl J., 1978,Mathematical Statistics (in Czech), SNTL/Alfa, Praha

3. Kunz J., Kovarik O., LauschmannH., Siegl J., Augustin P., 2010, Fractographic reconsti-
tution of fatigue crack growth in integrally stiffened panels,FATIGUE 2010, 2, 1711-1720

4. Kvasnička V., Posṕichal J., Tiňo P., 2000,Evolutionary Algorithms (in Slovak), STU Brati-
slava

5. Lauschmann H., Goldsmith N., 2009, Textural fractography of fatigue fractures, [In:] Fatigue
Crack Growth:Mechanics, Behavior and Prediction, A.F. Lignelli (Edit.), Nova SciencePublishers,
Inc., 125-166

6. Lauschmann H., Kunz J., Kovař́ik O., 2011, Analyzing crack growth rate in steel P92 (in
Czech), Research report V-KMAT-839/11, FNSPECTUPrague

7. Lauschmann H., Siegl J., Sumbera J., Siska F., Nedbal I., 2006, An unifying concept for
fatigue: The reference crack growth rate,Materials Characterization, 56, 257-265

8. Lauschmann H., Siska F., 2012, The reference texture: A proposal of a physical explanation,
International Journal of Fatigue, 43, 120-127

9. Mojzeš M., Klimt M., Kukal J., 2016, Feature selection via competitive levy flights, 2016
International Joint Conference on Neural Networks (IJCNN)

10. Mojzeš M., Kukal J., Lauschmann H., 2012, Sub-model testing in fractographic analy-
sis, Proceedings of Mendel 2012 Soft Computing Conference, Brno University Technology Press,
350-355

11. Nedbal I., Lauschmann H., Siegl J., Kunz J., 2008, Fractographic reconstitution of fatigue
crack history – Part II,Fatigue and Fracture of Engineering Materials and Structures, 31, 177-183

12. Ralston A., Rabinowitz P., 2001, A First Course in Numerical Analysis, Courier Dover Pu-
blications

13. SekeresovaZ., LauschmannH., 2008,Multi-fractal features of fatigue crack surfaces in relation
to crack growth rate,Materials Structure and Micromechanics of Fracture V, 567-568, 129-132

14. Schwarz G., 1978, Estimating the dimension of a model,Annals of Statistics, 6, 461-464

15. Tvrd́ik J., Křivý I., 2011, Hybrid adaptive differential evolution in partitional clustering, Pro-
ceedings of 17th International Conference on Soft Computing, Mendel 2011, VUT Brno, Brno,
1-8

16. Wilks S.S., 1962,Mathematical Statistics, JohnWiley & Sons, Inc., NewYork

17. Wooldridge J.M., 2002, Econometric Analysis of Cross Section and Panel Data, Cambridge,
MA:MIT Press

Manuscript received February 14, 2017; accepted for print May 16, 2017