title: Optimisation of matrix in R
body: <p>I'm new to optimisation/calibration of models in R, but i'm eager to learn and really need some help. My question relates to demographic modelling.</p>
  
  <p>I've done some research and found help <a href="https://stats.stackexchange.com/questions/136787/r-constrained-optimization-for-a-function-taking-a-matrix">here</a> and <a href="https://stackoverflow.com/questions/42421003/optimization-of-values-in-a-matrix">here</a> but neither have quite answered my question. </p>
  
  <p>I have a matrix of scalars (propensities) where each column must total to 1. These propensities are used to estimate the number of households that would arise from a given population (persons by age). The propensities model tends to overestimate the number of households in history (for which I know the true number of households). </p>
  
  <p>I want to calibrate the model to minimise the error in the number of households by tweaking the propensities such that the columns still add to 1 and propensities with an initial value of zero must remain zero.</p>
  
  <p>Simple example: </p>
  
  <pre><code>  # Propensities matrix
    mtx &lt;- matrix(c(0.00,  0.00,  0.85,  0.00,  0.15, 0.35,  0.45,  0.00,  
                 0.20,  0.00, 0.65,  0.15,  0.00,  0.20,  0.00), ncol = 3)
  
    # Population by age cohort
    pop &lt;- c(2600, 16200, 13400)
  
    # True number of households
    target &lt;- c(7000, 4500, 5500)
  
    # Function to optimise
    hh &lt;- function(mtx, pop, target) {
      # Estimate living arrangements
      x &lt;- mtx %*% pop 
      # Estimate number of households using parent cohorts (1,2 and 4)
      x &lt;- c(x[1,1]/2, x[2,1]/2, x[4,1]) - target
      return(x)
    }
  </code></pre>
  
  <p>I haven't included any of my code for the optimisation/calibration step as it would be embarrassing and I've haven't been able to get anything to work! </p>
  
  <p>Ideally i will have one set of propensities that generalises well for lots of different regions at the end of this process. Any advice on how i should go about achieving it? Helpful links?</p>
  
  <p><strong>Update</strong></p>
  
  <p>The snippet of code below executes the local search method as suggested by Enrico.</p>
  
  <pre><code>library(tidyverse)
  library(NMOF)
  
  data &lt;- list(mtx = matrix(c(0.00,  0.00,  0.90,  0.00,  0.10, 0.25,  0.50,  0.00,   
                              0.25,  0.00, 0.60,  0.20,  0.00,  0.20,  0.00), ncol = 3),
               pop = c(2600, 16200, 13400),
               target = c(7190, 4650, 5920))
  
  # True mtx
  mtx.true &lt;- matrix(c(0.00,  0.00,  0.75,  0.00,  0.25, 0.35,  0.45,  0.00,
                       0.20,  0.00, 0.65,  0.15,  0.00,  0.20,  0.00), ncol = 3)
  
  # Function to optimise
  households &lt;- function(x, data) {
  
    # Estimate living arrangements
    z &lt;- x %*% data$pop 
  
    # Estimate number of households using parent cohorts (1,2 and 4)
    z &lt;- c(z[1,1]/2, z[2,1]/2, z[4,1]) - data$target
    sum(abs(z))
  
  }
  
  # Local search function to perturb propensities
  neighbour &lt;- function(x, data) {
  
    # Choose random column from mtx
    i &lt;- sample(1:ncol(x), 1)
    # Select two non-zero propensities from mtx column
    j &lt;- which(x[, i] != 0) %&gt;% sample(2, replace = FALSE)
  
    # Randomnly select one to perturb positively 
    x[j[1], i] &lt;- 0.1 * (1 - x[j[1], i]) + x[j[1], i]
    # Perturb second propensity to ensure mtx column adds to 1
    x[j[2], i] &lt;- x[j[2], i] + (1 - sum(x[,i]))
  
    x
  
  }
  
  # Local search algorithm inputs 
  localsearch &lt;- list(x0 = data$mtx,             
                      neighbour = neighbour,
                      nS = 50000,                
                      printBar = FALSE)
  
  # Execute 
  now &lt;- Sys.time()
  solution &lt;- LSopt(OF = households, algo = localsearch, data)
  #&gt; 
  #&gt; Local Search.
  #&gt; Initial solution:  2695 
  #&gt; Finished.
  #&gt; Best solution overall: 425.25
  Sys.time() - now
  #&gt; Time difference of 6.33272 secs
  
  # Inspect propensity matrices
  print(solution$xbest)
  #&gt;           [,1]   [,2] [,3]
  #&gt; [1,] 0.0000000 0.3925  0.6
  #&gt; [2,] 0.0000000 0.4250  0.2
  #&gt; [3,] 0.2937976 0.0000  0.0
  #&gt; [4,] 0.0000000 0.1825  0.2
  #&gt; [5,] 0.7062024 0.0000  0.0
  print(mtx.true)
  #&gt;      [,1] [,2] [,3]
  #&gt; [1,] 0.00 0.35 0.65
  #&gt; [2,] 0.00 0.45 0.15
  #&gt; [3,] 0.75 0.00 0.00
  #&gt; [4,] 0.00 0.20 0.20
  #&gt; [5,] 0.25 0.00 0.00
  </code></pre>
  
  <p>Thanks!</p>