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?

  • 1
    $\begingroup$ 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? $\endgroup$
    – usεr11852
    Apr 15, 2018 at 12:28
  • $\begingroup$ 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? $\endgroup$
    – jo941
    Apr 15, 2018 at 12:40
  • $\begingroup$ Isn't this a convex problem? Why use simulated-annealing, which is designed for functions with many local minimum? $\endgroup$
    – shimao
    Apr 15, 2018 at 16:50
  • $\begingroup$ I guess, but my goal was to implement and compare several different optimization algorithms, just to get to know them better. $\endgroup$
    – jo941
    Apr 15, 2018 at 17:37

1 Answer 1


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))

            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)
[1] 530

[1] 1.018276

[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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.