tuning via maximizing likelihood

[Rein Halbersma]

Another important issue is "model mismatch", that is, when the true evaluation function (predicting the true w/l probabilities via the logistic function) and the heuristic evaluation function have different forms (e.g. E(x)=a+bx, H(c)=cx, where x is a scalar property of the position).

Both ML and LS will come up with a value for c, which will be different in both cases. Unfortunately there is no objective criterion (yet?) to decide which value is better. So this issue cannot be handled theoretically.

Model miss-specification aside, there is a theoretical sense in which MLE is superior to non-linear least squares (NLS).

For any finite-sample of N observations generated from a model with the "true" parameter vector, you can compute the vector of parameter estimates for both MLE and NLS. You can also do this for multiple of such finite samples, generating a sampling distribution of parameter estimates.

Under some mild conditions, both the MLE and NLS estimates converge in probability to the "true" parameter vector as the sample size N goes to infinity. In econometrics jargon, this means that both estimators are "consistent". In other words, the parameter estimates for MLE and NLS are in practice the same.

However, the variance of the sampling distribution of the parameter estimates scales at least as fast as 1/sqrt(N). Only for MLE is this bound strictly satisfied, and in econometrics jargon it is called the "most efficient" estimator of all consistent estimators. NLS always has a larger variance around the "true" parameter for finite samples.

Source: Micro-econometrics textbook by Cameron & Trivedi, chapter 5.

[Rein Halbersma]

LogL converges 3x faster (fewer iterations) than MSE if the positions have a win rate of 50% . For a win rate of 90% this advantage is ~9x, and for 99% win rate the advantage is ~35x. So the higher the win rate, the faster LogL converges compared to MSE. Note that in all cases the fitted parameter value is virtually identical.

While I also would like to believe that ML is superior to least squares I am not sure your test shows this. I think all you have been measuring is the performance of R's cg method for optimizing certain specific 1-variable functions.

I used the conjugate gradient method because in my application domain (draughts), the number of parameters is usually in the tens of thousands or even hundreds of thousands, which precludes methods that require a Hessian matrix of second-order derivatives (such as Newton's method). For similar reasons, gradient methods are the de facto standard in deep learning.

For the curious: the large number of parameters arises as follows. You can divide the 10x10 draughts board into 4 x 4 overlapping regions of 4x4 squares each (the overlap is 2 squares in each case). Each 4x4 region has 8 valid squares (draughts is played on the dark squares only) and can take on 3 valued (black, white or empty; kings are disregarded for positional patterns). Then one can compute an index in the range 1 ... 3^8 for each of the 16 regions.

In the eval, one simply computes the indices for each of the 16 patterns (depending on the board occupancy for the position being evaluated), and then one does 16 table lookups to get the corresponding parameter value. This is all fitted using either MSE or MLE, starting from initial values equal to zero. Some L1 or L2 regularization can also be used.

These type of evaluation functions were pioneered in draughts in Fabien Letouzey's draughts program Scan (3-times and reigning champion of the Computer Olympiad). Other patterns are of course also possible.

[Daniel Shawul]

We still need to determine A since log( L(Ax+b) ) = (1 - L(x)) * A, and this would mean using automatic differentiation via operator overloading or something like that. I was hoping there would be a surrogate minimizer of the ML objective similar to the one for the BT model, but i now see that this is impossible atleast without keeping track of factors for each parameter inside eval -- which i wanted to avoid.

[jdart]

Go programs also have gigantic numbers of evaluation parameters.

--Jon

[Daniel Shawul]

There was some typo in the sign of the draw score which I corrected below.

score is evaluation of position in centi-pawns relative to the side to move.
home_score is the score of being the side to move (not necessarily white)
draw_score is the threshold score of getting a draw

Rao-Kapper

Code: Select all

p(win:score) = logistic(score + home_score - draw_score) 
p(loss:score) = logistic(-score - home_score - draw_score) 
p(draw:score) = 1 - p(win:score) - p(loss:score) 

Glenn-David

Code: Select all

p(win:score) = gaussian(score + home_score - draw_score) 
p(loss:score) = gaussian(-score - home_score - draw_score) 
p(draw:score) = 1 - p(win:score) - p(loss:score) 

Davidson

Code: Select all

draw_gamma = score_to_gamma(draw_score)
f = draw_gamma * sqrt(logistic(score + home_score) *  logistic(-score - home_score) ) 
p(win:score) = logistic(score + home_score) / (1 + f) 
p(loss:score) = logistic(-score - home_score) / (1 + f) 
p(draw:score) = 1 - p(win:score) - p(loss:score) 

where:

Code: Select all

logistic(score) = 1 / (1 + pow(10.0,score / 400.0))
gaussian(score) = (1 + erf(-score / 400.0)) / 2
score_to_gamma(score) = pow(10.0,score / 400.0);

[Daniel Shawul]

Using a bayesian approach, we can add a prior distribution of parameters which acts like a regularization for the ML estimates. So we do a maximum a posteriori (MAP) estimate rather than maximum likelihood instead (ML). So far, I have not been able to get better results than my existing parameers for the eval terms, so it makes sense to keep that momentum ('or absence of it') going for a while by introducing a prior distribution over the evaluation parameters (theta)

Code: Select all

P(theta:results) = P(results:theta)P(theta)=likelihood*prior

[Rein Halbersma]

It would nice to understand this for evaluation functions depending on more parameters. The model would be a "true" (unknown) evaluation function E(x1,...,xn) predicting the correct w/l ratios (using the logistic function) depending on measurable parameters x1,..,xn and a heuristic evaluation function H(x1,...,xn).

I guess in first approximation E and H can be assumed to be linear in x1,...,xn. The challenge is to match up the coefficients of H with those of E.

In a realistic analysis draws should also be incorporated. The presence of draws may smooth out the effect of evaluation errors.

OK, I did some more thorough analyses. Below a code snippet in R (Works On My Machine).

What I do is generate a 4-term evaluation function. Two terms are 0/1 variables, generated by a binomial process. The other two are normally distributed variables.

I fix the "true" eval parameters to be -0.6190392 0.4054651 0.1000000 -0.1500000 respectively. I then generate linear predictors (called "eta", as in the literature on generalized models). I encode this as 0, 0.5 and 1 depending into which section between the theta thresholds the linear predictor falls. I also fix the "true" theta parameters to -1.386294 1.386294 (this corresponds to a 60% draw rate).

After that, I run two models: 1) Non-Linear Squares using a continuous LHS (y1 in my program) that has 0, 0.5 and 1 as values, and regress it on the eval features, 2) Cumulative Link Model using a one-hot encoded triple LHS (y3 in my program), and regress it on the eval features and also fitting the draw threshold parameters. I extract the fitted coefficients and compute the linear predictors (the evals used during search).

It turns out that there is perfect linear correlation between the two evals!! The NLS estimates are 65% smaller than the CLM estimates, in order to "squeeze" the eval score towards the value of a draw. The CLM estimates don't have to because this model already has the theta parameters as generalized intercepts.

Code: Select all

library(ordinal)
# make this reproducible 
set.seed(47110815) 

# Generate N logistically distributed values 
N = 1e5
x = cbind(rbinom(N, 1, .5), rbinom(N, 1, .5), rnorm(N), rnorm(N))
(w_star = c(qlogis(.35), qlogis(.6), .1, -.15))
theta_star = c(qlogis(.2), qlogis(.8))
eta = as.vector(x %*% w_star)
y_latent = eta + rlogis(N, 0, 1)

# 0/0.5/1 encoded
y1 = ifelse(y_latent > theta_star[2], 1, 
  ifelse(y_latent > theta_star[1], 1/2, 0)
)

# one-hot encoded of length 3
y3 = apply(cbind(y1 == 0, y1 == 1/2 , y1 == 1), c(1,2), as.integer)

# Use R library functions to get the most accurate results 
w0 = rep(0, length(w_star))

(nls1.est = nls(y1 ~ plogis(x %*% w), start = list(w = w0)))
(-logLik(nls1.est) / N) 
(clm1.est = clm(factor(y1) ~ x))
(-logLik(clm1.est) / N)

nls1.eta = x %*% coef(nls1.est)
clm1.eta = x %*% coef(clm1.est)[-(1:2)]

plot(clm1.eta, nls1.eta)
summary(lm(nls1.eta ~ clm1.eta))

Output:

Code: Select all

> (nls1.est = nls(y1 ~ plogis(x %*% w), start = list(w = w0)))
Nonlinear regression model
  model: y1 ~ plogis(x %*% w)
   data: parent.frame()
      w1       w2       w3       w4 
-0.39761  0.25896  0.06154 -0.10136 
 residual sum-of-squares: 9977

Number of iterations to convergence: 2 
Achieved convergence tolerance: 5.337e-06
> (-logLik(nls1.est) / N) 
'log Lik.' 0.2664951 (df=5)
> (clm1.est = clm(factor(y1) ~ x))
formula: factor(y1) ~ x

 link  threshold nobs  logLik    AIC       niter max.grad cond.H 
 logit flexible  1e+05 -94636.09 189284.17 6(0)  1.47e-12 1.1e+01

Coefficients:
      x1       x2       x3       x4 
-0.60832  0.40026  0.09369 -0.15488 

Threshold coefficients:
 0|0.5  0.5|1 
-1.380  1.387 
> (-logLik(clm1.est) / N)
'log Lik.' 0.9463609 (df=6)
> 
> nls1.eta = x %*% coef(nls1.est)
> clm1.eta = x %*% coef(clm1.est)[-(1:2)]
> 
> plot(clm1.eta, nls1.eta)
> summary(lm(nls1.eta ~ clm1.eta))

Call:
lm(formula = nls1.eta ~ clm1.eta)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0041136 -0.0010737  0.0000032  0.0010759  0.0034799 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
&#40;Intercept&#41; -1.467e-03  4.136e-06  -354.7   <2e-16 ***
clm1.eta     6.523e-01  9.861e-06 66150.5   <2e-16 ***
---
Signif. codes&#58;  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error&#58; 0.001266 on 99998 degrees of freedom
Multiple R-squared&#58;      1,	Adjusted R-squared&#58;      1 
F-statistic&#58; 4.376e+09 on 1 and 99998 DF,  p-value&#58; < 2.2e-16

For completeness (code not shown), I also reproduced the above R library functions with hand-written mean-squared error and log-likelihood formulas, and the outcomes agreed perfectly.

[Michel]

We still need to determine A since log( L(Ax+b) ) = (1 - L(x)) * A, and this would mean using automatic differentiation via operator overloading or something like that. I was hoping there would be a surrogate minimizer of the ML objective similar to the one for the BT model, but i now see that this is impossible atleast without keeping track of factors for each parameter inside eval -- which i wanted to avoid.

It may be more simple. If the evaluation function is linear in a parameter (are there good examples where it is not?) you can trivially differentiate it numerically with respect to that parameter. The result will be exact.

In this way one can even compute the Hessian of the objective function. Of course if you have a lot of parameters the Hessian will be a big matrix. But this is not such a big issue on modern computers. For example numpy has no problem whatsoever inverting a 1000x1000 matrix.

I guess if you need 36*3^8 parameters, as Rein was reporting for draughts, then using a second order method seems out of the question.

[Michel]

Using a bayesian approach, we can add a prior distribution of parameters which acts like a regularization for the ML estimates. So we do a maximum a posteriori (MAP) estimate rather than maximum likelihood instead (ML). So far, I have not been able to get better results than my existing parameers for the eval terms, so it makes sense to keep that momentum ('or absence of it') going for a while by introducing a prior distribution over the evaluation parameters (theta)

I do not think putting on a prior will make much of a difference, except perhaps for the speed of convergence.

What might be an interesting experiment is to make the draw_elo parameter (in the BE model) itself a linear combination of the features (i.e. it would vary from position to position). A model with a constant draw_elo parameter is obviously incorrect as in the endgame the expected draw rate should be much higher.

If implemented literally this would roughly double the number of parameters. There is an increased risk of over fitting however so one could cut this down by only including features that are expected to have an influence on the draw ratio, like game phase and perhaps king safety.

[jdart]

r * log( logistic(score) ) + (1 - r) log( 1 - logistic(score) )

This may have been covered, at least indirectly, but the issue I see with this function and draws is that the value of the loss is not zero when the predicted score is 0.5 and the actual score is 0.5. The loss function for draws does have its minimum value in that case. But it is not zero.

Each draw will increase the value of the objective even if it is predicted correctly, and this increase is more than what most loss and win positions contribute (because many will be more or less accurately predicted).

--Jon