# Simulated annealing for least squares/ linear regression

I am trying to optimize the solution of the linear regression problem by performing simulated annealing on the following loss function of least squares (MSE):

$$1/N (\sum (y_i-(ax_i+b))^2$$

I have studied a lot of examples of simulated annealing performed on functions with constant weights, but have not come across to something similar to that problem, which has a and b instead of constant numbers. Any ideas on how i could implement sa on that kind of equations? Should I somehow choose some constants to represent the a and b?

• Do you want to implement SAnn or just solve LR with SAnn? For example: myloss <- function(beta){sum((y-beta[1]-beta[2]*x )^2)}; N = 99; set.seed(N); x = runif(N); eps = rnorm(N); y = 1 + 3*x+ eps; optim(myloss, method= "SANN", par = c(2,2)) encodes a simple univariate LR task solve through SAnn. We can easily check this returns the same result as lm(y ~x), is that what you want? – usεr11852 Apr 15 '18 at 12:28
• I would like to find a global minimum of the function, by implementing Sann from scratch, to better familiarize myself with that method. In the code you wrote, beta is given as an argument, so I will just input random numbers as beta values? Or should I acquire the best solution somehow? – jo941 Apr 15 '18 at 12:40
• Isn't this a convex problem? Why use simulated-annealing, which is designed for functions with many local minimum? – shimao Apr 15 '18 at 16:50
• I guess, but my goal was to implement and compare several different optimization algorithms, just to get to know them better. – jo941 Apr 15 '18 at 17:37

Simulated annealing's fundamental mechanism is to compare the current state with a proposed state, then (probabilistically) decide whether to move to the proposed state or not. In your case, "state" = "parameter estimates", because those are what you are trying to find the optimum values for.

We can set up a state vector $\theta_i = (a^*_i,b^*_i)$, the estimates of the parameters at iteration $i$. You'll need to construct a proposal distribution which suggests a candidate $\theta_{i+1}$. One such proposal distribution could be a Student-t centered at the current values of the parameters with a fairly low shape parameter $\nu$, e.g., $\nu = 3$. (Having fat-tailed proposal distributions helps the algorithm get out of local minima.) The scale of the proposal distribution would have to be problem-specific.

The following code is a very simple implementation of SA for solving the univariate regression $y_i = a + bx_i + e_i$, and may help to get you started:

simulated_annealing <- function(theta, y, x, scale = 0.1, max_iter = 1000, step = 0.05) {
degrees_of_freedom <- length(y) - 2
temp <- 1

# Need to keep track of best state so far and current state
# We also keep track of the iteration count of the best state for tuning purposes
iter_best <- 1
theta_best <- theta
mse_best <- mse <- sum((y - theta[1] - theta[2]*x)^2) / degrees_of_freedom
for (k in 1:max_iter) {
temp <- temp * (1 - step)

# Generate proposals using a Student-t distribution
theta_prop <- theta + scale * rt(2, 3)
mse_prop <- sum((y - theta_prop[1] - theta_prop[2]*x)^2) / degrees_of_freedom

# Select the next state (either current or proposed)
if (mse_prop < mse || runif(1) < exp(-(mse_prop-mse)/temp)) {
theta <- theta_prop
mse <- mse_prop

if (mse_prop < mse_best) {
iter_best <- k
theta_best <- theta_prop
mse_best <- mse_prop
}
}
}
return(list(iter_best=iter_best, mse_best=mse_best, theta_best=theta_best))
}


Now to try it out, removing lots of extraneous text from the output:

> x <- rnorm(100)
> y <- 1 + 2*x + rnorm(100)
> # Compare to linear regression estimates
> summary(lm(y~x))

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.06306    0.10117   10.51   <2e-16 ***
x            2.07854    0.09896   21.00   <2e-16 ***

Residual standard error: 1.009 on 98 degrees of freedom
Multiple R-squared:  0.8182,    Adjusted R-squared:  0.8164

> # Starting values for SA parameter estimation: (0,0)
> simulated_annealing(c(0,0), y, x)
$iter_best [1] 530$mse_best
[1] 1.018276

\$theta_best
[1] 1.061634 2.083421


We observe that our parameter estimates are very, very close to those of the linear regression. Our RMSE is 1.009097, slightly higher than lm's 1.009082 that we can calculate by saving the lm results as an object and examining the contents.