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? 
 A: 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.
