Alternating optimization

The function call

The ao() function call with the default arguments looks as follows:

ao(
  f,
  initial,
  target = NULL,
  npar = NULL,
  gradient = NULL,
  ...,
  partition = "sequential",
  new_block_probability = 0.5,
  minimum_block_number = 2,
  minimize = TRUE,
  lower = -Inf,
  upper = Inf,
  iteration_limit = Inf,
  seconds_limit = Inf,
  tolerance_value = 1e-6,
  tolerance_parameter = 1e-6,
  tolerance_parameter_norm = function(x, y) sqrt(sum((x - y)^2)),
  base_optimizer = Optimizer$new("stats::optim", method = "L-BFGS-B"),
  verbose = FALSE
)

The arguments have the following meaning:

• f: The objective function to be optimized. By default, f is optimized over its first argument. If optimization should target a different argument or multiple arguments, use npar and target, see below. Additional arguments for f can be passed via the ... argument as usual.

• initial: Initial values for the parameters used in the AO procedure.

• gradient: Optional argument to specify the analytical gradient of f. If not provided, a finite-difference approximation will be used.

• partition: Specifies how parameters are partitioned for optimization. Can be one of the following:

  • "sequential": Optimizes each parameter block sequentially.

  • "random": Randomly partitions parameters in each iteration.

  • "none": No partitioning; equivalent to joint optimization.

  • Custom partition can be defined using a list of vectors of parameter indices.

• new_block_probability and minimum_block_number are only relevant if partition = "random". In this case, the former controls the probability for creating a new block, and the latter defines the minimum number of parameter blocks.

• minimize: Set to TRUE for minimization problems (default), or FALSE for maximization.

• lower and upper: Lower and upper limits for constrained optimization.

• iteration_limit is the maximum number of AO iterations before termination, while seconds_limit is the time limit in seconds. tolerance_value and tolerance_parameter (in combination with tolerance_parameter_norm) specify two other stopping criteria, namely when the change in the function value or the parameter distance, respectively, becomes smaller than these thresholds.

• base_optimizer: Numerical optimizer used for solving sub-problems.

• verbose: Set to TRUE to print status messages during the AO process.

A simple first example

The following is an implementation of the Himmelblau’s function \[f(x, y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2:\]

himmelblau <- function(x) (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2

This function has four identical local minima, for example in \(x = 3\) and \(y = 2\):

himmelblau(c(3, 2))
#> [1] 0

Minimizing Himmelblau’s function through alternating minimization for \(\mathbf{x}\) and \(\mathbf{y}\) with initial values \(\mathbf{x}^{(0)} = \mathbf{y}^{(0)} = 0\) can be accomplished as follows:

ao(f = himmelblau, initial = c(0, 0))
#> $estimate
#> [1]  3.584428 -1.848126
#> 
#> $value
#> [1] 9.606386e-12
#> 
#> $details
#>    iteration        value       p1        p2 b1 b2     seconds update_code
#> 1          0 1.700000e+02 0.000000  0.000000  0  0 0.000000000           0
#> 2          1 1.327270e+01 3.395691  0.000000  1  0 0.015664101           0
#> 3          1 1.743664e+00 3.395691 -1.803183  0  1 0.002752066           0
#> 4          2 2.847290e-02 3.581412 -1.803183  1  0 0.002337933           0
#> 5          2 4.687468e-04 3.581412 -1.847412  0  1 0.002980947           0
#> 6          3 7.368057e-06 3.584381 -1.847412  1  0 0.001838923           0
#> 7          3 1.164202e-07 3.584381 -1.848115  0  1 0.014685154           0
#> 8          4 1.893311e-09 3.584427 -1.848115  1  0 0.001490116           0
#> 9          4 9.153860e-11 3.584427 -1.848124  0  1 0.001165152           0
#> 10         5 6.347425e-11 3.584428 -1.848124  1  0 0.001085043           0
#> 11         5 9.606386e-12 3.584428 -1.848126  0  1 0.001051903           0
#> 
#> $seconds
#> [1] 0.04505134
#> 
#> $stopping_reason
#> [1] "change in function value is < 1e-06"

Here, we see the output of the alternating optimization procedure, which is a list that contains the following elements:

• estimate is the parameter vector at termination.

• value is the function value at termination.

• details is a data.frame with full information about the procedure: For each iteration (column iteration) it contains the function value (column value), parameter values (columns starting with p followed by the parameter index), the active parameter block (columns starting with b followed by the parameter index, where 1 stands for a parameter contained in the active parameter block and 0 if not), computation times in seconds (column seconds), and a code that summarizes whether the update got accepted (column update_code, where 0 stands for an accepted update, 1 for a rejected update due to an error when solving the sub-problem, and 2 for a rejected update because it did not improve the function value).

• seconds is the overall computation time in seconds.

• stopping_reason is a message why the procedure has terminated.

Using the analytical gradient

For the Himmelblau’s function, it is straightforward to define the analytical gradient as follows:

gradient <- function(x) {
  c(
    4 * x[1] * (x[1]^2 + x[2] - 11) + 2 * (x[1] + x[2]^2 - 7),
    2 * (x[1]^2 + x[2] - 11) + 4 * x[2] * (x[1] + x[2]^2 - 7)
  )
}

The gradient function will be used by ao() if defined through the gradient argument as follows:

ao(f = himmelblau, initial = c(0, 0), gradient = gradient)

The output is not shown here because it closely resembles the previous example, where the gradient was not specified and thus a finite-difference approximation was employed. However, in scenarios involving higher dimensions, utilizing the analytical gradient can notably improve both the speed and stability of the process.

Random parameter partitions

Another version of the AO procedure involves using a new, random partition of the parameters in every iteration. This approach can enhance the convergence rate and prevent being stuck in local optima. It is activated by setting partition = "random". The randomness can be adjusted using two parameters:

• new_block_probability determines the probability for creating a new block when building a new partition. Its value ranges from 0 (no blocks are created) to 1 (each parameter is a single block).

• minimum_block_number sets the minimum number of parameter blocks for random partitions. By default, it is configured to 2 to avoid generating trivial partitions.

The random partitions are build as follows:¹

procedure <- ao:::Procedure$new(
  npar = 10,
  partition = "random",
  new_block_probability = 0.5,
  minimum_block_number = 2
)
procedure$get_partition()
#> [[1]]
#> [1] 5
#> 
#> [[2]]
#> [1] 1 6 9
#> 
#> [[3]]
#> [1] 10
#> 
#> [[4]]
#> [1] 7
#> 
#> [[5]]
#> [1] 4 8
#> 
#> [[6]]
#> [1] 2 3

procedure$get_partition()
#> [[1]]
#> [1] 1 7 8
#> 
#> [[2]]
#> [1]  6 10
#> 
#> [[3]]
#> [1] 3 4
#> 
#> [[4]]
#> [1] 2
#> 
#> [[5]]
#> [1] 9
#> 
#> [[6]]
#> [1] 5

As an example of AO with random partitions, consider fitting a two-class Gaussian mixture model via maximizing the model’s log-likelihood function

\[\ell(\boldsymbol{\theta}) = \sum_{i=1}^n \log\Big( \lambda \phi_{\mu_1, \sigma_1^2}(x_i) + (1-\lambda)\phi_{\mu_2,\sigma_2^2} (x_i) \Big),\]

where the sum goes over all observations \(1, \dots, n\), \(\phi_{\mu_1, \sigma_1^2}\) and \(\phi_{\mu_2, \sigma_2^2}\) denote the normal density for the first and second cluster, respectively, and \(\lambda\) is the mixing proportion. The parameter vector to be estimated is thus \(\boldsymbol{\theta} = (\mu_1, \mu_2, \sigma_1, \sigma_2, \lambda)\). As there exists no closed-form solution for the maximum likelihood estimator \(\boldsymbol{\theta}^* = \arg\max_{\boldsymbol{\theta}} \ell(\boldsymbol{\theta})\), we need numerical optimization for finding the function optimum. The model is fitted to the following data:²

The following function calculates the log-likelihood value given the parameter vector theta and the observation vector data:³

normal_mixture_llk <- function(theta, data) {
  mu <- theta[1:2]
  sd <- exp(theta[3:4])
  lambda <- plogis(theta[5])
  c1 <- lambda * dnorm(data, mu[1], sd[1])
  c2 <- (1 - lambda) * dnorm(data, mu[2], sd[2])
  sum(log(c1 + c2))
}

The ao() call for performing alternating maximization with random partitions looks as follows, where we simplified the output for brevity:

out <- ao(
  f = normal_mixture_llk,
  initial = runif(5),
  data = datasets::faithful$eruptions,
  partition = "random",
  minimize = FALSE
)
round(out$details)
#>    iteration value p1 p2 p3 p4 p5 b1 b2 b3 b4 b5 seconds update_code
#> 1          0  -714  1  1  1  0  1  0  0  0  0  0       0           0
#> 2          1  -541  1  4  1  0  1  0  1  0  0  0       0           0
#> 3          1  -513  1  4  1  0  1  0  0  1  1  0       0           0
#> 4          1  -448  3  4  1  0  1  1  0  0  0  0       0           0
#> 5          1  -445  3  4  1  0  0  0  0  0  0  1       0           0
#> 6          2  -432  3  4  1  0  0  0  1  0  0  0       0           0
#> 7          2  -277  2  4 -1 -1 -1  1  0  1  1  1       0           0
#> 8          3  -276  2  4 -1 -1 -1  0  1  0  0  1       0           0
#> 9          3  -276  2  4 -1 -1 -1  1  0  1  1  0       0           0
#> 10         4  -276  2  4 -1 -1 -1  0  0  1  0  0       0           2
#> 11         4  -276  2  4 -1 -1 -1  1  1  0  1  1       0           0
#> 12         5  -276  2  4 -1 -1 -1  1  0  0  0  0       0           0
#> 13         5  -276  2  4 -1 -1 -1  0  1  1  1  1       0           0
#> 14         6  -276  2  4 -1 -1 -1  0  0  0  0  1       0           0
#> 15         6  -276  2  4 -1 -1 -1  1  1  1  1  0       0           0
#> 16         7  -276  2  4 -1 -1 -1  1  1  0  1  1       0           0
#> 17         7  -276  2  4 -1 -1 -1  0  0  1  0  0       0           0
#> 18         8  -276  2  4 -1 -1 -1  1  1  0  1  1       0           0
#> 19         8  -276  2  4 -1 -1 -1  0  0  1  0  0       0           0

More flexibility

The {ao} package offers some flexibility for performing AO.⁴

normal_mixture_llk <- function(data, mu, sd, lambda) {
  sd <- exp(sd)
  lambda <- plogis(lambda)
  c1 <- lambda * dnorm(data, mu[1], sd[1])
  c2 <- (1 - lambda) * dnorm(data, mu[2], sd[2])
  sum(log(c1 + c2))
}

ao(
  f = normal_mixture_llk,
  initial = runif(5),
  target = c("mu", "sd", "lambda"),
  npar = c(2, 2, 1),
  data = datasets::faithful$eruptions,
  partition = "random",
  minimize = FALSE
)

normal_mixture_llk <- function(mu, sd, lambda, data) {
  c1 <- lambda * dnorm(data, mu[1], sd[1])
  c2 <- (1 - lambda) * dnorm(data, mu[2], sd[2])
  sum(log(c1 + c2))
}
ao(
  f = normal_mixture_llk,
  initial = runif(5),
  target = c("mu", "sd", "lambda"),
  npar = c(2, 2, 1),
  data = datasets::faithful$eruptions,
  partition = "random",
  minimize = FALSE,
  lower = c(-Inf, -Inf, 0, 0, 0),
  upper = c(Inf, Inf, Inf, Inf, 1)
)