Optimización Monte Carlo
María Guzmán Martínez
2026
library(tinytex)
library(webshot2)
library(chromote)
library(plotly) # Gráfico dinámico
library(spFSR) # spsa_optimizer
rm(list=ls())Casella 1999
Introducción
De manera similar al problema de la integración, las diferencias entre el enfoque numérico y el enfoque de simulación para el problema \[ \max_{\theta \in \Theta} h(\theta) \tag{5.1.1} \] radican en el tratamiento de la función \(h\). (Note que (5.1.1) también abarca problemas de minimización al considerar \(-h\).) Al abordar un problema de optimización utilizando métodos numéricos deterministas, las propiedades analíticas de la función objetivo (convexidad, acotación, suavidad) suelen ser de primordial importancia. Para el enfoque de simulación, estamos más preocupados por \(h\) desde un punto de vista probabilístico (más que analítico). Obviamente, esta dicotomía es algo artificial, ya que existen enfoques de simulación donde no se utiliza la interpretación probabilística de \(h\). No obstante, el uso de las propiedades analíticas de \(h\) juega un papel menor en el enfoque de simulación.
Los métodos numéricos gozan de una historia más larga que los métodos de simulación (véase, por ejemplo, Kennedy y Gentle 1980 o Thisted 1988), pero los métodos de simulación han ganado atractivo debido a la relajación de las restricciones tanto en la regularidad del dominio como en la función \(h\). Por supuesto, puede existir un enfoque numérico alternativo que proporcione una solución exacta a (5.1.1), una propiedad rara vez alcanzada por un algoritmo estocástico, pero la simulación tiene la ventaja de omitir los pasos preliminares de diseñar un algoritmo y estudiar si se cumplen algunas condiciones de regularidad sobre \(h\). Esto es particularmente cierto cuando la función \(h\) es muy costosa de calcular.
Ejemplo 5.1.1 Procesamiento de señales. Ó Ruanaidh y Fitzgerald (1996) estudian datos de procesamiento de señales, de los cuales un modelo simple es (\(i = 1, \dots, N\)): \[ x_{i} = \alpha_{1}\cos(\omega t_{i}) + \alpha_{2}\sin(\omega t_{i}) + \epsilon_{i} \] \[ \epsilon_{i} \sim \mathcal{N}(0, \sigma^{2}), \] con parámetros desconocidos \(\alpha = (\alpha_{1}, \alpha_{2})\), \(\omega\), \(\sigma\), y tiempos de observación \(t_{1}, \dots, t_{N}\). La función de verosimilitud es entonces de la forma: \[ \sigma^{-N} \exp \left( -\frac{(x - G\alpha^{t})^{t}(x - G\alpha^{t})}{2\sigma^{2}} \right) \] con \(x = (x_{1}, \dots, x_{N})\) y \[ G = \begin{pmatrix} \cos(\omega t_{1}) & \sin(\omega t_{1}) \\ \vdots & \vdots \\ \cos(\omega t_{N}) & \sin(\omega t_{N}) \end{pmatrix} \] La distribución previa \(\pi(\alpha, \omega, \sigma) = \sigma^{-1}\) conduce a la distribución marginal: \[ \pi(\omega | x) \propto \left( x^{t}x - x^{t}G(G^{t}G)^{-1}G^{t}x \right)^{(2-N)/2} (\det G^{t}G)^{-1/2}, \tag{5.1.2} \] la cual, aunque explícita en \(\omega\), no es particularmente simple de calcular. Esta configuración también es ilustrativa de funciones con múltiples modas, como lo demuestran Ó Ruanaidh y Fitzgerald (1996).
Siguiendo a Geyer (1996), queremos distinguir entre dos enfoques de la optimización Monte Carlo:
Enfoque exploratorio: El objetivo es optimizar la función \(h\) describiendo todo su rango. Las propiedades reales de la función juegan un papel menor aquí, estando el aspecto Monte Carlo más estrechamente ligado a la exploración de todo el espacio \(\Theta\), aun cuando, por ejemplo, la pendiente de \(h\) pueda usarse para acelerar la exploración. (Tal técnica puede ser útil para describir funciones con múltiples modas, por ejemplo).
Enfoque basado en aproximación probabilística: Se basa en una aproximación probabilística de la función objetivo \(h\) y es en cierto modo un paso preliminar a la optimización real. Aquí, el aspecto Monte Carlo explota las propiedades probabilísticas de la función \(h\) para proponer una aproximación aceptable y se preocupa menos de explorar \(\Theta\). Veremos que este enfoque puede vincularse a métodos de datos faltantes, como el algoritmo EM.
Notamos también que Geyer (1996) solo considera que el segundo enfoque es “optimización Monte Carlo”. Obviamente, aunque estemos considerando estos dos enfoques diferentes por separado, podrían combinarse en un problema dado. De hecho, métodos como el algoritmo EM (5.3.3) o el algoritmo de Robbins-Monro (§5.5.3) aprovechan la aproximación Monte Carlo para mejorar su técnica de optimización particular.
Exploración Estocástica
Una Solución Básica
Hay una serie de casos donde el método de exploración es particularmente adecuado. En primer lugar, si \(\Theta\) está acotado, lo cual a veces se puede lograr mediante una reparametrización, un primer enfoque para la resolución de (5.1.1) es simular a partir de una distribución uniforme en \(\Theta\), \(u_{1}, \dots, u_{m} \sim \mathcal{U}_{\Theta}\), y usar la aproximación: \[ h_{m}^{*} = \max(h(u_{1}), \dots, h(u_{m})). \] Este método converge (cuando \(m\) tiende a \(\infty\)), pero puede ser muy lento ya que no tiene en cuenta ninguna característica específica de \(h\). Las distribuciones distintas de la uniforme, que posiblemente puedan estar relacionadas con \(h\), pueden funcionar mejor. En particular, en configuraciones donde la función de verosimilitud es extremadamente costosa de calcular, es mejor mantener al mínimo el número de evaluaciones de la función \(h\).
Esto conduce a una segunda dirección, a menudo más fructífera, que relaciona \(h\) con una distribución de probabilidad. Por ejemplo, si \(h\) es positiva y si: \[ \int_{\Theta} h(\theta) d\theta < +\infty, \] la resolución de (5.1.1) equivale a encontrar las modas de la densidad \(h\). Más generalmente, si no se cumplen estas condiciones, podemos transformar la función \(h(\theta)\) en otra función \(H(\theta)\) que cumpla con lo siguiente:
La función \(H\) es no negativa y satisface \(\int H < \infty\).
Las soluciones a (5.1.1) son aquellas que maximizan \(H(\theta)\) en \(\Theta\).
Por ejemplo, podemos tomar: \[ H(\theta) = \exp(h(\theta)/T) \quad \text{o} \quad H(\theta) = \frac{\exp\{h(\theta)/T\}}{1 + \exp\{h(\theta)/T\}} \] y elegir \(T\) para acelerar la convergencia o evitar máximos locales (como en el templado simulado; véase 5.2.3). Cuando el problema se expresa en términos estadísticos, resulta natural generar una muestra \((\theta_{1}, \dots, \theta_{m})\) a partir de \(h\) (o \(H\)) y aplicar un método estándar de estimación de la moda (o simplemente comparar los \(h(\theta_{i})\)). (En algunos casos, puede ser más útil descomponer \(h(\theta)\) como \(h(\theta) = h_{1}(\theta)h_{2}(\theta)\) y simular a partir de \(h_{1}\).)
El atractivo de la simulación es aún más claro en el caso de que \(h\) pueda representarse como: \[ h(\theta) = \int H(x,\theta) dx. \] En particular, si \(H(x,\theta)\) es una densidad y si es posible simular a partir de esta densidad, la solución de (5.1.1) es la moda de la distribución marginal de \(\theta\). (Aunque esta configuración pueda parecer forzada o incluso artificial, veremos que incluye el caso de modelos de datos faltantes.)
A continuación, examinamos varios métodos para encontrar máximos que pueden clasificarse como métodos exploratorios.
Métodos del Gradiente
Como se mencionó en la 1.4, el método del gradiente es un enfoque numérico determinista para el problema (5.1.1). Produce una secuencia \((\theta_{j})\) que converge a la solución exacta de (5.1.1), \(\theta^{*}\), cuando el dominio \(\Theta \subset \mathbb{R}^{d}\) y la función \((-h)\) son ambos convexos. La secuencia \((\theta_{j})\) se construye de manera recursiva mediante: \[ \theta_{j+1} = \theta_{j} + \alpha_{j} \nabla h(\theta_{j}), \quad \alpha_{j} > 0, \tag{5.2.1} \] donde \(\nabla h\) es el gradiente de \(h\). Para varias elecciones de la secuencia \((\alpha_{j})\) (véase Thisted 1988), el algoritmo converge al máximo (único).
En configuraciones más generales (es decir, cuando la función o el espacio son menos regulares), la ecuación (5.2.1) se puede modificar mediante perturbaciones estocásticas para lograr de nuevo la convergencia, como se describe en detalle en Rubinstein (1981) o Duflo (1996, pp. 61-63). Una de estas modificaciones estocásticas consiste en elegir una segunda secuencia \((\beta_{j})\) para definir la cadena \((\theta_{j})\) mediante: \[ \theta_{j+1} = \theta_{j} + \frac{\alpha_{j}}{2\beta_{j}} \Delta h(\theta_{j}, \beta_{j}\zeta_{j}) \zeta_{j}. \tag{5.2.2} \] Las variables \(\zeta_{j}\) se distribuyen uniformemente en la esfera unitaria \(\|\zeta\| = 1\) y: \[ \Delta h(x,y) = h(x+y) - h(x-y) \] aproxima \[2\|y\| \nabla h(x)\]. A diferencia del enfoque determinista, este método no avanza necesariamente a lo largo de la pendiente más pronunciada en \(\theta_{j}\), pero esta propiedad a veces es una ventaja en el sentido de que puede evitar quedar atrapado en máximos locales o puntos de silla de \(h\).
La convergencia de \((\theta_{j})\) a la solución \(\theta^{*}\) depende nuevamente de la elección de \((\alpha_{j})\) y \((\beta_{j})\). Notamos de paso que \((\theta_{j})\) se puede ver como una cadena de Markov no homogénea (véase la Definición 4.2.3) que converge casi con certeza a un valor dado. El estudio de estas cadenas es bastante complicado debido a su núcleo de transición cambiante (véase Winkler 1995 para algunos resultados en esta dirección). Sin embargo, condiciones suficientemente fuertes como la disminución de \(\alpha_{j}\) hacia 0 y de \(\alpha_{j}/\beta_{j}\) hacia una constante no nula son suficientes para garantizar la convergencia de la secuencia \((\theta_{j})\).
Calculo de máximos
Ejemplo 1
En R, existen varias funciones para resolver problemas de optimización. La más simple es la función optimize, que procesa objetivos unidimensionales.
Ejemplo 5.1. Al considerar la maximización de la verosimilitud de una muestra de una distribución de Cauchy \(C(\theta, 1)\),
\[\ell(\theta \vert x_1, \dots, x_n) = \prod_{i=1}^n \frac{1}{1 + (x_i - \theta)^2}\]
la secuencia de máximos (es decir, de los estimadores de máxima verosimilitud o MLEs) converge a \(\theta^* = 0\) cuando \(n\) tiende a \(\infty\).
xm=rcauchy(500, location = 0, scale = 1)
f=function(theta){-sum(log(1+(x-theta)^2))}
for (i in 1:500) {
x = xm[1:i]
mi = optimise(f, interval = c(-10, 10), maximum = TRUE)$maximum
}
mi## [1] 0.0004780747donde la log-verosimilitud se maximiza secuencialmente a medida que
la muestra aumenta. Sin embargo, al mirar directamente a la
verosimilitud, el uso de optimise eventualmente produce una
secuencia divergente debido a que la verosimilitud se vuelve demasiado
pequeña alrededor de las \(n = 300\)
observaciones, a pesar de que las secuencias de MLEs son idénticas hasta
este punto para ambas funciones.
Ejemplo 2
El método del gradiente es un enfoque numérico determinista para optimización. En el contexto de encontrar un máximo, la regla de actualización iterativa se define como: \[ \theta_{j+1} = \theta_j + \alpha_j \nabla h(\theta_j), \quad \alpha_j > 0 \] donde \(\nabla h(\theta)\) representa el gradiente de la función objetivo \(h\) y \(\alpha_j\) es la tasa de aprendizaje.
Consideremos la función cóncava de dos variables orientada a maximizar: \[ h(\theta_1, \theta_2) = -(\theta_1 - 2)^2 - (\theta_2 - 3)^2 + 10 \] El punto óptimo (máximo global) de esta función se encuentra analíticamente en \(\theta^* = (2, 3)\), donde el valor máximo es \(h(2,3) = 10\).
Cálculo del Gradiente: Calculo de las derivadas parciales de primer orden respecto a cada parámetro: \[ \nabla h(\theta) = \begin{bmatrix} \frac{\partial h}{\partial \theta_1} \\ \frac{\partial h}{\partial \theta_2} \end{bmatrix} = \begin{bmatrix} -2(\theta_1 - 2) \\ -2(\theta_2 - 3) \end{bmatrix} \]
# función a maximizar: h
h <- function(theta) {return(-((theta[1] - 2)^2) - ((theta[2] - 3)^2) + 10)}# Delta h
grad_h <- function(theta) {
dt1 <- -2 * (theta[1] - 2)
dt2 <- -2 * (theta[2] - 3)
return(c(dt1, dt2))
}# Algoritmo del Gradiente (Ascenso)
ascenso_gradiente <- function(theta_inicial, alpha = 0.1, tol = 1e-6, max_iter = 1000) {
theta <- theta_inicial
for (i in 1:max_iter) {
grad <- grad_h(theta)
# Actualización de parámetros
theta_nuevo <- theta + alpha * grad
# Criterio de convergencia
if (sqrt(sum((theta_nuevo - theta)^2)) < tol) {
cat("Convergencia alcanzada en la iteración:", i, "\n")
break
}
theta <- theta_nuevo
}
return(theta)
}
# Punto de inicio
theta_inicio <- c(0, 0)
tasa_aprendizaje <- 0.1
resultado <- ascenso_gradiente(theta_inicio, alpha = tasa_aprendizaje)## Convergencia alcanzada en la iteración: 62resultado # Punto óptimo encontrado (theta)## [1] 1.999998 2.999996 h(resultado) # Valor máximo de la función## [1] 10Gráfico del resultado:
Función para encontrar el máximo usando optim
# función a maximizar: h
h <- function(theta) {
theta1 <- theta[1]
theta2 <- theta[2]
return(-((theta1 - 2)^2) - ((theta2 - 3)^2) + 10)
}
# optim por defecto minimiza, luego con fnscale = -1 se maximiza
resultado <- optim(
par = c(0, 0), # cóndiciones iniciales
fn = h, # Función objetivo
method = "BFGS", # Método basado en gradiente (Quasi-Newton)
control = list(fnscale = -1) # Invierte el sentido para buscar un MÁXIMO
)
resultado$par # Punto óptimo encontrado (theta)## [1] 2 3resultado$value # Valor máximo de la función## [1] 10Nota: El algoritmo BFGS (Método Quasi-Newton) es uno de los métodos Quasi-Newton más eficientes para resolver problemas de optimización no lineal sin restricciones. A diferencia del método de Newton puro, que requiere calcular la matriz Hessiana exacta (\(\nabla^2 f(x)\)) en cada iteración —lo cual es computacionalmente costoso—, el método BFGS aproxima la inversa de la matriz Hessiana de forma iterativa utilizando únicamente evaluaciones de la función y sus gradientes.
Fórmula de actualización de la matriz inversa: En cada iteración \(k\), se mantiene una matriz de aproximación \(B_k \approx (\nabla^2 f(x_k))^{-1}\). La fórmula matemática para actualizar esta matriz del paso \(k\) al paso \(k+1\) es:
\[ B_{k+1} = (I - \rho_k s_k y_k^T) B_k (I - \rho_k y_k s_k^T) + \rho_k s_k s_k^T \]
Donde:
\(s_k = x_{k+1} - x_k\): Vector de cambio en la posición de los parámetros.
\(y_k = \nabla f(x_{k+1})-\nabla f(x_k)\): Vector de cambio en el gradiente.
\(\rho_k = \dfrac{1}{y_k^T s_k}\): Factor de escala escalar.
\(I\): Matriz identidad de dimensiones concordantes.
Pasos del Algoritmo
-Inicialización: Se selecciona un punto inicial \(x_0\) y se inicializa la matriz aproximada \(B_0\) comúnmente como la matriz identidad (\(I\)).
Dirección de búsqueda (\(p_k\)): Se determina la dirección en la que se moverá el algoritmo utilizando la matriz actual y el gradiente: \[ p_k = -B_k \nabla f(x_k) \]
Búsqueda lineal (Line Search): Se calcula una longitud de paso óptima \(\alpha_k\) en la dirección \(p_k\) para garantizar el descenso (o ascenso) adecuado de la función.
Actualización de la posición: \[ x_{k+1} = x_k + \alpha_k p_k \]
Actualización de la matriz: Se calculan \(s_k\) e \(y_k\) con los nuevos valores y se aplica la fórmula BFGS para obtener \(B_{k+1}\).
Criterio de parada: El proceso se repite hasta que la norma del gradiente sea menor a una tolerancia establecida (\(\|\nabla f(x_k)\| < \text{tol}\)) o se alcance el número máximo de iteraciones.
Ejemplo 3
Del ejercicion, 5.2.2: Se tiene \(\mathbb{R}^{2}\):
\[ h(x,y) = (x \sin(20y) + y \sin(20x))^{2} \cosh(\sin(10x)x) + (x \cos(10y) - y \sin(10x))^{2} \cosh(\cos(20y)y), \]
h=function(x,y){(x*sin(20*y)+y*sin(20*x))^2*cosh(sin(10*x)
*x)+(x*cos(10*y)-y*sin(10*x))^2*cosh(cos(20*y)*y)}
# 1. Definir la malla de puntos para x e y
# Reducimos un poco el rango si deseas ver el detalle central con más claridad,
# o puedes ajustarlo según lo necesites (ej. de -2 a 2)
x <- seq(-2, 2, length.out = 150)
y <- seq(-2, 2, length.out = 150)
# 2. Crear la función h(x, y)
h <- function(x, y) {
# Como outer() trabaja con matrices, usamos funciones vectorizadas de R
term1 <- (outer(x, rep(1, length(y)), function(xi, yi) xi * sin(20 * yi)) +
outer(rep(1, length(x)), y, function(xi, yi) yi * sin(20 * xi)))^2 *
cosh(outer(x, rep(1, length(y)), function(xi, yi) sin(10 * xi) * xi))
term2 <- (outer(x, rep(1, length(y)), function(xi, yi) xi * cos(10 * yi)) -
outer(rep(1, length(x)), y, function(xi, yi) yi * sin(10 * xi)))^2 *
cosh(outer(rep(1, length(x)), y, function(xi, yi) cos(20 * yi) * yi))
return(term1 + term2)
}
# 3. Evaluar la función en la malla 2D
z <- outer(x, y, h)
# (Opcional) Limitar valores extremadamente altos de z si el Cosh crece demasiado rápido,
# para evitar que la gráfica se distorsione visualmente:
z[z > 50] <- NA
# 4. Generar la gráfica 3D interactiva con Plotly
plot_ly(x = x, y = y, z = z) |>
add_surface(
colorscale = "Viridis", # Puedes cambiar a "Plasma", "Jet", "Cividis", etc.
showscale = TRUE
) |>
layout(
title = "Gráfica de h(x, y)",
scene = list(
xaxis = list(title = "x"),
yaxis = list(title = "y"),
zaxis = list(title = "h(x, y)"),
camera = list(
eye = list(x = 1.5, y = 1.5, z = 1.2)
)
)
)Encontrar el máximo de la función con el algoritmo BFGS:
# función h
h <- function(v)
{
x <- v[1]
y <- v[2]
val <- (x*sin(20*y)+y*sin(20*x))^2*cosh(sin(10*x)
*x)+(x*cos(10*y)-y*sin(10*x))^2*cosh(cos(20*y)*y)
return(val)
}
set.seed(200)
# optim con el método BFGS
resultado <- optim(
par = c(0.65, 0.8), # Punto inicial de theta_0, evitar el [0, 0] porque el gradiente puede anularse
fn = h,
method = "BFGS",
control = list(fnscale = -1) # fnscale = -1 para maximizar
)
# Resultados
resultado$par # máximos## [1] 445.1801 1964.6822resultado$value # Valor máximo de h## [1] 1.180784e+153Ejemplo 4
Ejemplo 5.2.: La verosimilitud asociada con el modelo de mezcla
\[(5.2) \quad \frac{1}{4} \mathcal{N}(\mu_1, 1) + \frac{3}{4} \mathcal{N}(\mu_2, 1)\]
es bimodal. Calcular los máximos de la función.
ATENCIÓN: Las funciones multimodales tienen varios picos (máximos locales); luego se le debe dar un valor inicial (par) cercano al pico que deseas encontrar.
mu1 <- -2 # Media 1
mu2 <- 3 # Media 2
sigma1 <- 1 # Desviación estándar 1
sigma2 <- 1 # Desviación estándar 2
h <- function(x) {
0.25 * dnorm(x, mean = mu1, sd = sigma1) +
0.75 * dnorm(x, mean = mu2, sd = sigma2)
}
resultado_pico1 <- optim(
par = 1, # Valor inicial cerca del primer pico
fn = h,
method = "BFGS",
control = list(fnscale = -1) # fnscale = -1 para maximizar
)
resultado_pico1$par## [1] 2.999994resultado_pico1$value # h evaluada en el punto## [1] 0.2992071Cálculo de mínimos
La función nlm de R que busca el mínimo de una función basándose en el método de Newton-Raphson:
\[\theta_{i+1} = \theta_i - \left[ \frac{\partial^2 h}{\partial \theta \partial \theta^T} (\theta_i) \right]^{-1} \frac{\partial h}{\partial \theta} (\theta_i)\]
donde la matriz de las segundas derivadas se llama Hessiana y el vector de las primeras derivadas es el gradiente (a veces denotado como \(\nabla h\)). Este método es perfecto cuando \(h\) es cuadrática, pero puede deteriorarse cuando \(h\) es altamente no lineal, y obviamente no funciona cuando el dominio \(\Theta\) es irregular. También depende claramente del punto de partida \(\theta_0\) cuando \(h\) tiene varios mínimos.
Ejemplo 1
Cálculo de mínimo de una función (pag. 250 del libro): \[ f(x,y)=x^2+y^2 \]
# Definir la malla de puntos para x e y
x <- seq(-3, 3, length.out = 50)
y <- seq(-3, 3, length.out = 50)
# Crear la función f(x, y) = x^2 + y^2
f <- function(x, y) {
x^2 + y^2
}
# Evaluar la función en la malla
z <- outer(x, y, f)
# Generar la gráfica 3D interactiva
plot_ly(x = x, y = y, z = z) |>
add_surface() |>
layout(
title = "Gráfica de f(x,y) = x^2 + y^2",
scene = c(
xaxis = list(title = "X"),
yaxis = list(title = "Y"),
zaxis = list(title = "Z")
)
)De acuerdo con Marsden, el mínimo está en el punto \((0,0)\)
# función
h<- function(theta) {
theta[1]^2 + theta[1]^2
}
# optim por defecto minimiza
resultado <- optim(
par = c(3, 3), # Punto inicial de búsqueda
fn = h, # Función objetivo
method = "BFGS", # Método basado en gradiente (Quasi-Newton)
control = list(fnscale = 1) # Encuentra el mínimo
)
resultado$par # Punto mínimo## [1] 5.492764e-16 3.000000e+00resultado$value # valor de la función en el punto encontrado## [1] 1.508406e-29El método BFGS se queda atrapado en un mínimo local
Ejemplo 6
Minimización de h en \(\mathbb{R}^{2}\): \[ h(x,y) = (x \sin(20y) + y \sin(20x))^{2} \cosh(\sin(10x)x) + (x \cos(10y) - y \sin(10x))^{2} \cosh(\cos(20y)y), \]
cuyo mínimo global es \(0\), alcanzado en los valores \((x,y) = (0,0)\). Dado que esta función tiene muchos mínimos locales, no satisface las condiciones bajo las cuales se garantiza que los métodos de minimización estándar proporcionen el mínimo global.
Por otro lado, la distribución en \(\mathbb{R}^{2}\) con densidad proporcional a \(\exp(-h(x,y))\) se puede simular, aunque no sea una distribución estándar, utilizando, por ejemplo, técnicas de Monte Carlo basadas en cadenas de Markov, y se puede derivar una aproximación convergente del mínimo de \(h(x,y)\) a partir del mínimo de los \(h(x_{i},y_{i})\) resultantes. Una alternativa es simular a partir de la densidad proporcional a:
\[ h_{1}(x,y) = \exp \{ -(x \sin(20y) + y \sin(20x))^{2} - (x \cos(10y) - y \sin(10x))^{2} \}, \]
lo que elimina el cálculo tanto de \(\cosh\) como de \(\sinh\) en el paso de simulación.
La exploración puede ser particularmente difícil cuando el espacio no es conexo (o tal vez ni siquiera conexo). En tales casos, la simulación de una muestra \((\theta_{1}, \dots, \theta_{m})\) puede ser mucho más rápida que un método numérico aplicado a (5.1.1).
Podemos aplicar la construcción iterativa (5.2.2) a la función multimodal \(h(x,y)\) con diferentes secuencias de \(\alpha_{j}\) y \(\beta_{j}\). La Figura 5.2.2 y la Tabla 5.2.1 ilustran que, dependiendo del valor inicial, el algoritmo converge a diferentes mínimos locales de la función \(h\). Aunque hay ocasiones en las que la secuencia \(h(\theta_{j})\) aumenta y evita algunos mínimos locales, las soluciones son bastante distintas para las tres secuencias diferentes, tanto en ubicación como en valores.
El número de iteraciones necesarias para lograr la estabilidad de \(\theta_{T}\) también varía con la elección de \((\alpha_{j}, \beta_{j})\). Tenga en cuenta que el Caso 1 da como resultado una evaluación muy deficiente del mínimo, ya que la rápida disminución de \((\alpha_{j})\) se asocia con grandes saltos en las primeras iteraciones. El Caso 2 converge a los mínimos locales más cercanos, y el Caso 3 ilustra una característica general del método del gradiente estocástico: tasas de disminución más lentas de la secuencia \((\alpha_{j})\) tienden a alcanzar mejores mínimos. También es destacable la convergencia final a lo largo de un valle de \(h\) después de algunos saltos grandes iniciales.
Este enfoque todavía está bastante cercano a los métodos numéricos en el sentido de que requiere un conocimiento preciso de la función \(h\), que no necesariamente está disponible.
Trayectorias del gradiente estocástico para tres elecciones diferentes de las secuencias \((\alpha_{j})\) y \((\beta_{j})\) y punto inicial \((0.65, 0.8)\) para la misma secuencia \((\zeta_{j})\) en (5.2.2). Los niveles de gris son tales que las tonalidades más oscuras significan elevaciones más altas. La función \(h\) a minimizar se define en el Ejemplo 5.2.1.
- \(\alpha_{j} = \beta_{j} = 1/(10j)\)
- \(\alpha_{j} = \beta_{j} = 1/(100j)\)
- \(\alpha_{j} = 1/(10\log(1+j))\), \(\beta_{j} = 1/j\)
Resultados de tres ejecuciones de gradiente estocástico para la minimización de la función \(h\) del Ejemplo 5.2.1 con diferentes valores de \((\alpha_{j}, \beta_{j})\) y punto de inicio \((0.65, 0.8)\). La iteración \(T\) se obtiene mediante la regla de parada \(\|\theta_{T} - \theta_{T-1}\| < 10^{-5}\).
| \(\alpha_{j}\) | \(\beta_{j}\) | \(\theta_{T}\) | \(h(\theta_{T})\) | \(\min_{t} h(\theta_{t})\) | Iteración \(T\) |
|---|---|---|---|---|---|
| \(1/(10j)\) | \(1/(10j)\) | \((-0.166, 1.02)\) | \(1.287\) | \(0.115\) | \(50\) |
| \(1/(100j)\) | \(1/(100j)\) | \((0.629, 0.786)\) | \(0.00013\) | \(0.00013\) | \(93\) |
| \(1/(10\log(1+j))\) | \(1/j\) | \((0.0004, 0.245)\) | \(4.24 \times 10^{-6}\) | \(2.163 \times 10^{-7}\) | \(58\) |
Código para encontrar un mínimo de h:
set.seed(200)
# minimizar la función h
h <- function(v) {
x <- v[1]
y <- v[2]
val <- (x*sin(20*y)+y*sin(20*x))^2*cosh(sin(10*x)
*x)+(x*cos(10*y)-y*sin(10*x))^2*cosh(cos(20*y)*y)
}
resultado <- optim(
par = c(0.65, 0.8), # Punto de inicio theta_0, evitar el [0, 0] porque el gradiente puede anularse
fn = h,
method = "BFGS",
control = list(maxit = 50, # Número de iteraciones
fnscale = 1 # Encuentra el mínimo
)
)
# Resultados
resultado$par # Coordenadas óptimas encontradas## [1] 0.6283185 0.7853981resultado$value # mínimo global/local de h## [1] 1.86981e-12Código para encontrar el mínimo de h:
# 1. Definición de la función objetivo h(x,y)
h <- function(theta) {
x <- theta[1]
y <- theta[2]
term1 <- (x * sin(20 * y) + y * sin(20 * x))^2 * cosh(sin(10 * x) * x)
term2 <- (x * cos(10 * y) - y * sin(10 * x))^2 * cosh(cos(20 * y) * y)
return(term1 + term2)
}
# 2. Generar UNA ÚNICA secuencia de ruidos zeta_j para garantizar
# que todos los casos compartan la misma trayectoria aleatoria base (Figura 5.2.2)
set.seed(1234) # Semilla ajustada para direcciones estables
max_steps <- 1000
phi_seq <- runif(max_steps, 0, 2 * pi)
zeta_matrix <- cbind(cos(phi_seq), sin(phi_seq))
# 3. Función de simulación usando la secuencia fija de zeta
stochastic_gradient_fixed <- function(case, theta0 = c(0.65, 0.8), tol = 1e-5) {
theta <- theta0
history_theta <- matrix(NA, nrow = max_steps, ncol = 2)
history_h <- numeric(max_steps)
history_theta[1, ] <- theta
history_h[1] <- h(theta)
j <- 1
converged <- FALSE
while (j < max_steps && !converged) {
if (case == 1) {
alpha_j <- 1 / (10 * j)
beta_j <- 1 / (10 * j)
} else if (case == 2) {
alpha_j <- 1 / (100 * j)
beta_j <- 1 / (100 * j)
} else if (case == 3) {
alpha_j <- 1 / (10 * log(1 + j))
beta_j <- 1 / j
}
# Usar el zeta de la secuencia compartida
zeta_j <- zeta_matrix[j, ]
# Diferencias finitas para estimar el gradiente
h_plus <- h(theta + beta_j * zeta_j)
h_minus <- h(theta - beta_j * zeta_j)
delta_h <- h_plus - h_minus
# Dirección de descenso de gradiente (minimización)
step <- (alpha_j / (2 * beta_j)) * delta_h * zeta_j
# Evitar saltos desproporcionados fuera de la región [-2, 2]
step_norm <- sqrt(sum(step^2))
if (step_norm > 0.5) {
step <- (step / step_norm) * 0.5
}
theta_next <- theta - step
j <- j + 1
history_theta[j, ] <- theta_next
history_h[j] <- h(theta_next)
# Criterio de parada: ||theta_T - theta_{T-1}|| < 10^-5
if (sqrt(sum((theta_next - theta)^2)) < tol) {
converged <- TRUE
}
theta <- theta_next
}
T_iter <- j
return(list(
theta_T = history_theta[T_iter, ],
h_theta_T = history_h[T_iter],
min_h_t = min(history_h[1:T_iter]),
T_iter = T_iter
))
}
# 4. Ejecución de los 3 casos
res1 <- stochastic_gradient_fixed(1)
res2 <- stochastic_gradient_fixed(2)
res3 <- stochastic_gradient_fixed(3)
# 5. Imprimir tabla con formato similar al libro
results_table <- data.frame(
Caso = c(1, 2, 3),
theta_T = c(
sprintf("(%.3f, %.3f)", res1$theta_T[1], res1$theta_T[2]),
sprintf("(%.3f, %.3f)", res2$theta_T[1], res2$theta_T[2]),
sprintf("(%.4f, %.3f)", res3$theta_T[1], res3$theta_T[2])
),
`h(theta_T)` = c(
sprintf("%.3f", res1$h_theta_T),
sprintf("%.5f", res2$h_theta_T),
sprintf("%.2e", res3$h_theta_T)
),
`min_h(theta_t)` = c(
sprintf("%.3f", res1$min_h_t),
sprintf("%.5f", res2$min_h_t),
sprintf("%.2e", res3$min_h_t)
),
Iteracion_T = c(res1$T_iter, res2$T_iter, res3$T_iter)
)
print(results_table)## Caso theta_T h.theta_T. min_h.theta_t. Iteracion_T
## 1 1 (0.314, 0.470) 0.000 0.000 32
## 2 2 (0.628, 0.785) 0.00000 0.00000 13
## 3 3 (0.0000, 0.180) 6.02e-09 6.02e-09 50Funciones del R para encontrar máximos: noisySPSA y spFSR
Función: spFSR
# Función del Ejemplo 5.2.1
h <- function(theta) {
x <- theta[1];
y <- theta[2]
(x * sin(20 * y) + y * sin(20 * x))^2 * cosh(sin(10 * x) * x) +
(x * cos(10 * y) - y * sin(10 * x))^2 * cosh(cos(20 * y) * y)
}# Función genérica SPSA (Simultaneous Perturbation Stochastic Approximation)
spsa_optimizer <- function(f, theta0, n_iter = 100, a = 0.1, c = 0.1, alpha = 0.602, gamma = 0.101) {
p <- length(theta0)
theta <- theta0
# Constante A recomendada por Spall (10% de las iteraciones totales)
A <- 0.1 * n_iter
history_theta <- matrix(NA, nrow = n_iter, ncol = p)
history_f <- numeric(n_iter)
for (k in 1:n_iter) {
# Secuencias de decaimiento a_k y c_k
a_k <- a / (k + A)^alpha
c_k <- c / k^gamma
# Perturbación simétrica (+1 o -1 con probabilidad 0.5)
delta <- sample(c(-1, 1), size = p, replace = TRUE)
# Evaluación de la función en dos puntos perturbados
theta_plus <- theta + c_k * delta
theta_minus <- theta - c_k * delta
f_plus <- f(theta_plus)
f_minus <- f(theta_minus)
# Estimación del gradiente estocástico
ghat <- (f_plus - f_minus) / (2 * c_k * delta)
# Actualización de theta (resta para minimizar)
theta <- theta - a_k * ghat
history_theta[k, ] <- theta
history_f[k] <- f(theta)
}
return(list(
par = theta,
value = f(theta),
history_theta = history_theta,
history_f = history_f
))
}
# Valores encontrados
res <- spsa_optimizer(f = h, theta0 = c(0.65, 0.8), n_iter = 200)
res$par # Mínimo encontrado (theta)## [1] 0.6549647 0.4914398res$value # Valor de h(theta)## [1] 6.617323e-05Método del recocido
Cálculo de máximos
Ejemplo 1
Dada la función \[ h(x,y) = (x \sin(20y) + y \sin(20x))^{2} \cosh(\sin(10x)x) + (x \cos(10y) - y \sin(10x))^{2} \cosh(\cos(20y)y) \] encontrar el máximo.
Algoritmo SANN
# función h
h <- function(v) {
x <- v[1]
y <- v[2]
(x*sin(20*y)+y*sin(20*x))^2*cosh(sin(10*x)
*x)+(x*cos(10*y)-y*sin(10*x))^2*cosh(cos(20*y)*y)
}
set.seed(200)
# Optimización estocástica
res_sann <- optim(
par = c(0.65, 0.8), # Punto de inicio
fn = h, # Función objetivo a maximizar
method = "SANN", # Simulated Annealing: método del Recocido Simulado
control = list(fnscale = -1)
)
res_sann$par # punto máximo (theta_T)## [1] -1.104258 5.027524res_sann$value # Valor de h(theta_T)## [1] 2860.416Ejemplo 2
Función: \[ f(x, y, z) = x^2 + y^2 + z^2 \]
# 1. Definir una función para generar los puntos de una esfera de radio 'r'
crear_esfera <- function(r, color_hex, nombre_nivel) {
# Crear malla de ángulos (theta y phi)
theta <- seq(0, 2 * pi, length.out = 40) # Longitud
phi <- seq(0, pi, length.out = 40) # Latitud
# Coordenadas cartesianas
x <- outer(cos(theta), sin(phi), "*") * r
y <- outer(sin(theta), sin(phi), "*") * r
z <- matrix(cos(phi), nrow = length(theta), ncol = length(phi), byrow = TRUE) * r
list(x = x, y = y, z = z, color = color_hex, name = nombre_nivel)
}
# 2. Definir los valores de las superficies de nivel (c = 1, 2, 3 por ejemplo)
# Nota: r es la raíz cuadrada de c (ya que x^2 + y^2 + z^2 = r^2 = c)
niveles <- list(
list(c = 1, r = 1, color = "rgba(31, 119, 180, 0.4)"),
list(c = 2, r = sqrt(2), color = "rgba(255, 127, 14, 0.4)"),
list(c = 3, r = sqrt(3), color = "rgba(44, 160, 44, 0.4)")
)
# 3. Construir la gráfica interactiva con Plotly
p <- plot_ly()
for (nivel in niveles) {
esfera <- crear_esfera(nivel$r, nivel$color, paste("c =", nivel$c))
p <- add_surface(p,
x = esfera$x,
y = esfera$y,
z = esfera$z,
surfacecolor = matrix(nivel$r, 40, 40),
colorscale = list(c(0, 1), c(nivel$color, nivel$color)),
name = paste("c =", nivel$c),
showscale = FALSE)
}
# 4. Configurar el diseño de los ejes y el título
p <- p |> layout(
title = "Superficies de nivel de f(x,y,z) = x^2 + y^2 + z^2",
scene = list(
xaxis = list(title = "X"),
yaxis = list(title = "Y"),
zaxis = list(title = "Z")
)
)
# Mostrar la gráfica
pEjemplo 4
Función: \[ f(x, y, z) = x^2 + y^2 - z^2 \]
# 1. Función para generar las mallas del hiperboloide de una hoja
crear_hiperboloide <- function(c_val, color_hex, nombre_nivel) {
# Rango para el radio r (debe empezar en sqrt(c_val) si c_val > 0, o en 0 si c_val <= 0)
r_min <- ifelse(c_val > 0, sqrt(c_val), 0)
r <- seq(r_min, r_min + 3, length.out = 30)
theta <- seq(0, 2 * pi, length.out = 30)
grid <- expand.grid(r = r, theta = theta)
# Coordenadas (la rama positiva y negativa de z)
# Usamos sqrt(grid$r^2 - c_val)
z_vals <- sqrt(pmax(0, grid$r^2 - c_val))
# Generar tanto la parte superior como la inferior del hiperboloide
x_sup <- grid$r * cos(grid$theta)
y_sup <- grid$r * sin(grid$theta)
z_sup <- z_vals
x_inf <- grid$r * cos(grid$theta)
y_inf <- grid$r * sin(grid$theta)
z_inf <- -z_vals
list(
list(x = matrix(x_sup, 30, 30), y = matrix(y_sup, 30, 30), z = matrix(z_sup, 30, 30), color = color_hex, name = paste(nombre_nivel, "(Sup)")),
list(x = matrix(x_inf, 30, 30), y = matrix(y_inf, 30, 30), z = matrix(z_inf, 30, 30), color = color_hex, name = paste(nombre_nivel, "(Inf)"))
)
}
# 2. Definir algunos valores de c para las superficies de nivel (ej: c = -2, -1, 0, 1, 2)
valores_c <- c(-2, -1, 0, 1, 2)
colores <- c("rgba(31, 119, 180, 0.3)",
"rgba(255, 127, 14, 0.3)",
"rgba(44, 160, 44, 0.3)",
"rgba(214, 39, 40, 0.3)",
"rgba(148, 103, 189, 0.3)")
# 3. Construir la gráfica interactiva con Plotly
p <- plot_ly()
for (i in seq_along(valores_c)) {
c_val <- valores_c[i]
hojas <- crear_hiperboloide(c_val, colores[i], paste("c =", c_val))
# Añadir la parte superior e inferior de cada nivel
for (hoja in hojas) {
p <- add_surface(p,
x = hoja$x,
y = hoja$y,
z = hoja$z,
colorscale = list(c(0, 1), c(hoja$color, hoja$color)),
showscale = FALSE,
name = hoja$name)
}
}
# 4. Configurar el diseño de los ejes y el título
p <- p |> layout(
title = "Superficies de nivel de f(x,y,z) = x^2 + y^2 - z^2",
scene = list(
xaxis = list(title = "X"),
yaxis = list(title = "Y"),
zaxis = list(title = "Z")
)
)
# Mostrar la gráfica
pMás gráficos, pag 90 del pdf.
Cálculo de mínimos
Ejemplo 1
Función:
\[ f(x,y)=x^2+y^2 \] Calcular el mínimo global.
h <- function(theta) {
theta[1]^2 + theta[2]^2
}
# Optimización global usando Recocido Simulado (SANN)
resultado_global <- optim(
par = c(10, 10), # Punto inicial lejano para probar la búsqueda global
fn = h, # Función objetivo
method = "SANN", # Método estocástico para búsqueda global
control = list(
maxit = 20000, # Número de iteraciones para explorar mejor
fnscale = 1 # fnscale = 1 para minimizar
)
)
resultado_global$par # Coordenadas del punto mínimo global## [1] -0.004128356 0.004275348resultado_global$value # Valor mínimo global ## [1] 3.532192e-05Ejemplo 2
\[ f(x,y)=x^2-y^2 \]
h <- function(theta) {
theta[1]^2 - theta[2]^2
}
# Optimización global usando Recocido Simulado (SANN)
resultado_global <- optim(
par = c(5, 5), # Punto inicial lejano para probar la búsqueda global
fn = h, # Función objetivo
method = "SANN", # Método estocástico para búsqueda global
control = list(
maxit = 20000, # Número de iteraciones para explorar mejor
fnscale = 1 # fnscale = 1 para minimizar
)
)
resultado_global$par # Coordenadas del punto mínimo global## [1] -8.097942 901.690567resultado_global$value # Valor mínimo global ## [1] -812980.3Ejemplo 3
Función: \[ f(x,y) = x^{1/3} y^{1/3} \]
x <- seq(-10, 10, length.out = 100)
y <- seq(-10, 10, length.out = 100)
# 3. Crear una función para la raíz cúbica real
raiz_cubica <- function(t) {
sign(t) * abs(t)^(1/3)
}
# 4. Calcular los valores de z sobre la cuadrícula (matriz)
z <- outer(x, y, function(x, y) raiz_cubica(x) * raiz_cubica(y))
# 5. Generar el gráfico interactivo en 3D
plot_ly(x = ~x, y = ~y, z = ~z) %>%
add_surface(
colorscale = "Viridis",
contours = list(
z = list(show = TRUE, usecolormap = TRUE, highlightcolor = "#ff0000", project = list(z = TRUE))
)
) %>%
layout(
title = "Gráfica de f(x,y) = x^(1/3) * y^(1/3)",
scene = list(
xaxis = list(title = "Eje X"),
yaxis = list(title = "Eje Y"),
zaxis = list(title = "Eje Z")
)
)Ejemplo 4
\[ z = \cos(xy), \quad -3 \leq x \leq 3, \, -3 \leq y \leq 3 \]
# Definir intervalos para x e y
x_b <- seq(-3, 3, length.out = 100)
y_b <- seq(-3, 3, length.out = 100)
# Calcular matriz z
z_b <- outer(x_b, y_b, function(x, y) cos(x * y))
# Graficar
plot_ly(x = ~x_b, y = ~y_b, z = ~z_b) %>%
add_surface(colorscale = "Viridis") %>%
layout(
title = "Gráfica (b): z = cos(xy)",
scene = list(xaxis = list(title = "X"), yaxis = list(title = "Y"), zaxis = list(title = "Z"))
)Ejemplo 6
\[ z = \frac{(x^2 - y^2)}{(x^2 + y^2)}, \quad 0 < |x| \leq 1.0, \, |y| \leq 1 \]
# Definir intervalos excluyendo el centro exacto para evitar NaN
x_c <- seq(-1, 1, length.out = 100)
y_c <- seq(-1, 1, length.out = 100)
# Calcular matriz z
z_c <- outer(x_c, y_c, function(x, y) {
num <- x^2 - y^2
den <- x^2 + y^2
res <- num / den
res[x == 0 & y == 0] <- NA # Evitar indeterminación en el origen
return(res)
})
# Graficar
plot_ly(x = ~x_c, y = ~y_c, z = ~z_c) %>%
add_surface(colorscale = "Plasma") %>%
layout(
title = "Gráfica (c): z = (x^2 - y^2) / (x^2 + y^2)",
scene = list(xaxis = list(title = "X"), yaxis = list(title = "Y"), zaxis = list(title = "Z"))
)Ejemplo 8
\[ z = \frac{[\operatorname{sen}(2x^2 + 3y^2)]}{(x^2 + y^2)}, \quad 0 < |x| \leq 3.0, \, |y| \leq 3 \]
# Definir intervalos para x e y
x_d <- seq(-3, 3, length.out = 150)
y_d <- seq(-3, 3, length.out = 150)
# Calcular matriz z evitando la división por cero en (0,0)
z_d <- outer(x_d, y_d, function(x, y) {
num <- sin(2 * x^2 + 3 * y^2)
den <- x^2 + y^2
res <- num / den
res[x == 0 & y == 0] <- NA # Evitar indeterminación
return(res)
})
# Graficar superficie en 3D
plot_ly(x = ~x_d, y = ~y_d, z = ~z_d) %>%
add_surface(colorscale = "Viridis") %>%
layout(
title = "Gráfica (d): z = sin(2x^2 + 3y^2) / (x^2 + y^2)",
scene = list(xaxis = list(title = "X"), yaxis = list(title = "Y"), zaxis = list(title = "Z"))
)Ejemplo 9
$$ z = , < |x| , , |y|
$$
# Definir intervalos para x e y
x_e <- seq(-3, 3, length.out = 150)
y_e <- seq(-3, 3, length.out = 150)
# Calcular matriz z evitando la división por cero en (0,0)
z_e <- outer(x_e, y_e, function(x, y) {
num <- sin(x * y)
den <- x^2 + y^2
res <- num / den
res[x == 0 & y == 0] <- NA # Evitar indeterminación
return(res)
})
# Graficar superficie en 3D
plot_ly(x = ~x_e, y = ~y_e, z = ~z_e) %>%
add_surface(colorscale = "Plasma") %>%
layout(
title = "Gráfica (e): z = sin(xy) / (x^2 + y^2)",
scene = list(xaxis = list(title = "X"), yaxis = list(title = "Y"), zaxis = list(title = "Z"))
)Ejemplo 10
Dada la distribución pareto, encontrar sus mínimos.
Notas
Considerar el caso especial \(f: U \subset \mathbb{R}^n \to \mathbb{R}\). Aquí \(Df(x)\) es una matriz de \(1 \times n\): \[Df(x) = \left[ \frac{\partial f}{\partial x_1} \dots \frac{\partial f}{\partial x_n} \right]\] Formamos el correspondiente vector \((\partial f / \partial x_1, \dots, \partial f / \partial x_n)\), llamado el gradiente de \(f\) y denotado por \(\operatorname{grad} f\) o \(\nabla f\).
Si \(f: U \subset \mathbb{R}^3 \to \mathbb{R}\) es diferenciable, el gradiente de \(f\) en \((x, y, z)\) es el vector en el espacio \(\mathbb{R}^3\) dado por: \[\operatorname{grad} f = \left( \frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z} \right)\]
Este vector también se denota por \(\nabla f\) o \(\nabla f(x, y, z)\). Así, \(\nabla f\) es simplemente la matriz de las derivadas \(Df\), escrita como vector.
TEOREMA: Suponer que \(\nabla f(\mathbf{x}) \neq \mathbf{0}\). Entonces \(\nabla f(\mathbf{x})\) apunta en la dirección a lo largo de la cual \(f\) crece más rápido.
Si queremos movernos en una dirección en la cual \(f\) va a crecer más rápidamente, debemos proceder en la dirección \(\nabla f(x)\). Análogamente, si deseamos movernos en una dirección en la cual \(f\) decrece más rápido, deberemos proceder en la dirección \(-\nabla f(x)\).
Sea \(f: U \subset \mathbb{R}^n \to \mathbb{R}\) una función escalar dada.
Un punto \(x_0 \in U\) se llama mínimo local de \(f\) si existe una vecindad \(V\) de \(x_0\) tal que para todos los puntos \(x \in V\), \[f(x) \ge f(x_0)\].
\(x_0 \in U\) es un máximo local si existe una vecindad \(V\) de \(x_0\) tal que \[f(x) \le f(x_0)\] para todo \(x \in V\).
El punto \(x_0 \in U\) es un extremo local o relativo, si es mínimo local o máximo local.
Un punto \(x_0\) es un punto crítico de \(f\) si \[Df(x_0) = 0\]. Un punto crítico que no es un extremo local se llama punto silla.
Teorema: Si \(U \subset \mathbb{R}^n\) es abierto, \(f: U \subset \mathbb{R}^n \to \mathbb{R}\) es diferenciable y \(\mathbf{x}_0 \in U\) es un extremo local, entonces \(Df(\mathbf{x}_0) = 0\); esto es, \(\mathbf{x}_0\) es un punto crítico de \(f\).
DEFINICIÓN: Sea \(A \subset \mathbb{R}^n\). Un punto \(\mathbf{x} \in \mathbb{R}^n\) es punto frontera de \(A\) si toda vecindad de \(\mathbf{x}\) contiene al menos un punto en \(A\) y al menos un punto fuera de \(A\).
DEFINICIÓN: Suponer que \(f: U \subset \mathbb{R}^n \to \mathbb{R}\) tiene derivadas parciales de segundo orden \((\partial^2 f / \partial x_i \partial x_j)(\mathbf{x}_0)\), para \(i, j = 1, \dots, n\), en un punto \(\mathbf{x}_0 \in U\). El hessiano de \(f\) en \(\mathbf{x}_0\) es la función cuadrática definida por \[ Hf(\mathbf{x}_0)(h) = \frac{1}{2} \sum_{i,j=1}^{n} \frac{\partial^2 f}{\partial x_i \partial x_j}(\mathbf{x}_0) h_i h_j. \]
TEOREMA: Si \(f: U \subset \mathbb{R}^n \to \mathbb{R}\) es de clase \(C^2\), \(\mathbf{x}_0 \in U\) es un punto crítico de \(f\) y el hessiano \(Hf(\mathbf{x}_0)\) es definitivamente positivo, entonces \(\mathbf{x}_0\) es un mínimo relativo de \(f\). De manera análoga, si \(Hf(\mathbf{x}_0)\) es definitivamente negativo, entonces \(\mathbf{x}_0\) es un máximo relativo.
En realidad, probaremos que los extremos son estrictos. Se dice que un máximo relativo \(\mathbf{x}_0\) es estricto si \(f(\mathbf{x}) < f(\mathbf{x}_0)\) para \(\mathbf{x}\) cercano, con \(\mathbf{x} \neq \mathbf{x}_0\). De manera análoga se define un mínimo relativo estricto.
DEFINICIÓN: Suponer que \(f: A \to \mathbb{R}\) es una función definida en un conjunto \(A\) de \(\mathbb{R}^2\) o \(\mathbb{R}^3\). Se dice que un punto \(\mathbf{x}_0 \in A\) es un punto de máximo absoluto (o de mínimo absoluto) de \(f\) si \(f(\mathbf{x}) \le f(\mathbf{x}_0)\) (o \(f(\mathbf{x}) \ge f(\mathbf{x}_0)\)) para todo \(\mathbf{x} \in A\).
DEFINICIÓN: Se dice que un conjunto \(D \subset \mathbb{R}^n\) es acotado si existe un número \(M > 0\) tal que \(\|\mathbf{x}\| < M\) para todo \(\mathbf{x} \in D\). Un conjunto es cerrado si contiene a todos sus puntos frontera.
TEOREMA DEL MÁXIMO Y DEL MÍNIMO: Sea \(D\) cerrado y acotado en \(\mathbb{R}^n\) y sea \(f: D \to \mathbb{R}\) continua. Entonces \(f\) alcanza sus valores máximo y mínimo en algunos puntos \(\mathbf{x}_0\) y \(\mathbf{x}_1\) de \(D\).