Regularización estadística de problemas inversos incorrectos para ellos. Turchin (parte 1)

Hola Habr! Hoy queremos contarte qué hace el laboratorio de métodos de experimentos de física nuclear , parte de JetBrains Research .

¿Dónde están JetBrains y dónde está la física nuclear? Acordamos sobre la base del amor por Kotlin, aunque en esta publicación no hablaremos de él. Nuestro grupo se enfoca en el desarrollo de software de análisis, modelado y escritura de datos para científicos, y por lo tanto se enfoca en la cooperación y el intercambio de conocimientos con las empresas de TI.

En este artículo queremos hablar sobre el método de regularización estadística que popularizamos , propuesto por V.F. Turchin en los años 70 del siglo XX, y su implementación en forma de código en Python y Julia.

La presentación será bastante detallada, por lo que aquellos que tienen claros los problemas inversos pueden ir directamente a los ejemplos y leer la teoría en este artículo .


Ocurrencia del problema: ¿por qué alguien debería regularizarse?


Si es suficiente ignorar, cualquier medición en el experimento se puede describir de la siguiente manera: hay un cierto dispositivo que captura el espectro o la señal de un proceso y muestra algunos números de acuerdo con los resultados de la medición. Nuestra tarea como investigadores, al observar estos números y conocer la estructura del dispositivo, es comprender cuál era el espectro o la señal medidos. Es decir, frente a lo que se llama el problema inverso . Si imagina esto matemáticamente, obtenemos esta ecuación (que, por cierto, se llama la ecuación de Fredholm del primer tipo ):

f(y)=abdxK(x,y)φ(x)

De hecho, esta ecuación describe lo siguiente: nuestro dispositivo de medición está representado aquí por su función de hardware K(x,y), que actúa sobre el espectro estudiado u otra señal de entradaφ , por el cual el investigador observa la señal de salidaf(y) . El objetivo del investigador es restaurar la señalφ por conocidof(y) yK(x,y) . También puede formular esta expresión en forma de matriz, reemplazando funciones con vectores y matrices:

fm=Kmnφn

Parecería que la reconstrucción de señales no es una tarea difícil, ya que tanto la ecuación de Fredholm como el sistema de ecuaciones lineales (incluso sobredeterminadas) tienen una solución exacta. Entonces probémoslo. Deje que la señal medida se describa como la suma de dos gausses:

φ(x)=2N(2,0.16)+N(4,0.04)

Como instrumento, tomamos el integrador más simple: una matriz que traduce nuestra señal en una suma acumulativa utilizando la función Heaviside:

Kmn=θ(xmyn)

El tipo de la señal medida y nuestro dispositivo, así como el resultado de la medición se muestran en el gráfico.

Es importante que cualquier medición real tenga un error, por lo que estropearemos un poco nuestro resultado agregando ruido normal, lo que da un error de medición del cinco por ciento.

Restauraremos la señal por el método de mínimos cuadrados:

φ=(KTK)1KTf

Y como resultado obtenemos:

En realidad, sobre esto, podríamos terminar el artículo, una vez más convencidos de la impotencia de los métodos idealistas de las matemáticas frente a la dura y dura realidad física, e ir a fumar soldadores.

Pero primero, averigüemos qué causó que nos ocurriera esta falla. Obviamente, el punto son los errores de medición, pero ¿a qué afectan? El hecho es que Jacques Hadamard (el mismo que agregó un guión a la fórmula Cauchy - Hadamard) dividió todas las tareas en correctas e incorrectas.

Recordando los clásicos: “No tiene sentido buscar una solución, si la hay. Se trata de cómo lidiar con una tarea que no tiene solución. Esta es una pregunta profundamente fundamental ... "- no hablaremos sobre las tareas correctas e inmediatamente tomaremos las incorrectas. Afortunadamente, ya hemos cumplido esto: la ecuación de Fredholm escrita anteriormente es un problema inverso incorrecto, incluso con fluctuaciones infinitesimales en los datos de entrada (e incluso nuestros errores de medición están lejos de ser infinitesimales), la solución de la ecuación obtenida de manera analítica exacta puede diferir arbitrariamente de la verdadera .

Puede leer la prueba de esta declaración en el primer capítulo del trabajo clásico del académico A.N. Tikhonova "Métodos para resolver problemas mal planteados". Este libro tiene consejos sobre qué hacer con tareas incorrectas, pero la técnica descrita allí tiene una serie de inconvenientes que se corrigen en el método Turchin. Pero primero, describimos los principios generales de trabajar con tareas incorrectas: ¿qué hacer si se encuentra con esa tarea?

Dado que la tarea en sí no puede ofrecernos nada, tendremos que cometer un pequeño delito: completar la tarea con datos para que sea correcta, en otras palabras, ingrese alguna * información a priori adicional sobre la tarea * (este proceso se llama regularización de la tarea). En contraste con el método clásico de Tikhonov, basado en la introducción de funciones de regularización parametrizadas, el método Turchin recurre, por otro lado, a los métodos bayesianos.

Descripción teórica de la regularización estadística


Estrategia


Formulamos nuestro problema en términos de estadística matemática: de acuerdo con una implementación bien conocida f(que medimos en el experimento) necesitamos estimar el valor del parámetroφ . FuncionalS^ computaciónφ basadof llamaremosestrategia. Para determinar qué estrategias son más óptimas, introducimos unafunción de pérdida cuadrática. La función de pérdida real puede ser cualquiera, ¿por qué elegimos la cuadrática? Debido a que cualquier función de pérdida cercana a su mínimo puede ser aproximada por una función cuadrática:

L(φ,S^[f])=||φ^S^[f])||L2,

Dónde - la mejor solución. Luego, las pérdidas para nuestra estrategia elegida están determinadaspor la función de riesgo:φ^

RS^[f](φ)E[L(φ,S^[f])]=L(φ,S^[f])P(f|φ)df.

aquí P(f|φ)determina la densidad de probabilidad de nuestro conjunto sobre el cual se promedia la pérdida. Este conjunto está formado por una hipotética repetición múltiple de mediciones.f para un determinadoφ . De este modo,P(f|φ) es la densidad de probabilidad que conocemosf obtenido en el experimento.

Según el enfoque bayesiano, se propone considerarφ como unavariable aleatoriacondensidad de probabilidad a priori P(φ), expresando lafiabilidad devarias soluciones a nuestro problema.P(φ) se determina en base a la información que existe antes del experimento. Entonces, la elección de la estrategia óptima se basa en minimizar elriesgo a posteriori:

rS^(φ)EφEf[L(φ,S^[f])|φ]

En este caso, la estrategia óptima es bien conocida:

S^[f]=E[φ|f]=φP(φ|f)dφ,

donde es la densidad posterior P(φ|f)está determinado por el teorema de Bayes:

P(φ|f)=P(φ)P(f|φ)dφP(φ)P(f|φ)

Este enfoque nos permitirá determinar la varianza (función de correlación) de la solución resultante:

D(x1,x2)=E[φ(x1)S^[f](x1)][φ(x2)S^[f](x2)]

Entonces, hemos obtenido la solución óptima a nuestro problema al introducir una densidad a priori P(φ). ¿Podemos decir algo sobre ese mundo de funciones?φ(x)que viene dada por la densidad a priori?

Si la respuesta a esta pregunta es no, entonces tendremos que aceptar todos los posiblesφ(x)igualmente probable y volver a una solución irregular. Por lo tanto, debemos responder a esta pregunta afirmativamente.

En esto consiste precisamente el método de regularización estadística: la regularización de una solución mediante la introducción de información adicional a priori sobreφ(x). Si el investigador ya tiene información a priori (densidad a prioriP(φ)), solo puede calcular la integral y obtener la respuesta.

Si no existe dicha información, el siguiente párrafo describe qué información mínima puede tener un investigador y cómo usarla para obtener una solución regularizada.

Información a priori


Como han demostrado los científicos británicos, en el resto del mundo les gusta diferenciarse. Además, si el matemático hace preguntas sobre la legalidad de esta operación, el físico cree de manera optimista que las leyes de la naturaleza se describen mediante funciones "buenas", es decir, suaves.

En otras palabras, lo hace más suaveφ(x)mayor densidad de probabilidad a priori. Así que tratemos de introducir una probabilidad a priori basada en la suavidad. Para hacer esto, recordamos que la introducción de información a priori es algo de violencia contra el mundo, lo que obliga a las leyes de la naturaleza a verse de una manera conveniente para nosotros.

Esta violencia debe minimizarse, y al introducir una densidad de probabilidad a priori, es necesario que la información de Shannon sobreφ(x)contenida en P(φ)fue mínimo Formalizando lo anterior, derivamos la forma de una densidad a priori basada en la suavidad de la función. Para hacer esto, buscaremos un extremo de información condicional:

I[P(φ)]=lnP(φ)P(φ)dφmin

Bajo las siguientes condiciones:

  1. Condición de suavidad φ(x). PermitirΩEs una matriz determinada que caracteriza la suavidad de la función. Luego requerimos que se logre un cierto valor de la suavidad funcional:

    (φ,Ωφ)P(φ)dφ=ω


    ω. .
  2. : P(φ)dφ=1
    :

    Pα(φ)=αRg(Ω)/2detΩ1/2(2π)N/2exp(12(φ,αΩφ))

    αω, , , , . α, . :

    1. α
    2. () α, α
    3. αP(α|f). Este enfoque es correcto porque proporciona una buena aproximación de la integral si los datos experimentales contienen suficiente información sobreα.

El primer caso es de poco interés para nosotros. En el segundo caso, debemos calcular una integral tan fea aquí:

φi=dφφiP(f|φ)dαP(α)αRg(Ω)2exp(α2(φ,Ωφ))dφP(f|φ)dαP(α)αRg(Ω)2exp(α2(φ,Ωφ))

Para el tercer caso, podemos obtener el valor de la integral analíticamente para los ruidos gaussianos en el experimento (esto se considerará a través de la sección).

También se debe tener en cuenta que no hemos usado en ningún lado, queΩEs un operador de suavidad. De hecho, aquí podemos usar cualquier otro operador (o su combinación lineal), solo la suavidad de la función es la forma más obvia de información a priori que podemos usar.

Muestreo


Hablamos de funciones, pero cualquier dispositivo real no puede medir no solo un continuo, sino incluso un conjunto contable de puntos. Siempre tomamos medidas en un conjunto finito de puntos, por lo tanto, nos vemos obligados a llevar a cabo procedimientos de discretización y transición de la ecuación integral a la matriz. En el método de regularización estadística, procedemos de la siguiente manera: descompondremosφ(x) sobre algún sistema de funciones {Tn}:

φ(x)=nφnTn(x).

Por lo tanto, los coeficientes de esta expansión forman algún vector φque es un vector en el espacio funcional.

Como espacio funcional, podemos tomar el espacio de Hilbert o, por ejemplo, el espacio de polinomios. Además, la elección de la base en estos espacios está limitada solo por su imaginación (intentamos trabajar con las series trigonométricas de Fourier, poligandra y splines cúbicas).

Entonces los elementos de la matrizK calculado como:

Kmn=(K^Tn(x))(ym),

Dónde ym- puntos en los que se realizaron mediciones. Elementos matricialesΩ calcularemos por la fórmula:

Ωij=ab(dpTi(x)dx)(dpTj(x)dx)dx,

Dónde ay b- los límites del intervalo en el que se define la función φ(x).

Para volver a calcular los errores, use la fórmula de dispersión de una combinación lineal de variables aleatorias:

D[φ(x)]=D[nφnTn(x)]=i,jφiφjcov(Ti(x),Tj(x)).

Debe tenerse en cuenta que, en algunos casos, la representación de una función utilizando un vector de dimensión finita conduce a una pérdida parcial o cambio de información. De hecho, podemos considerar la algebraización como un tipo de regularización, aunque débil e insuficiente para convertir una tarea incorrecta en una correcta. Pero, de todos modos, ahora hemos pasado de buscarφ(x) para buscar vectores φy en la siguiente sección lo encontramos.

Estuche de ruido gaussiano


El caso en que los errores en el experimento se distribuyen de acuerdo con Gauss es notable
porque se puede obtener una solución analítica a nuestro problema. Dado que la información y los errores a priori tienen una forma gaussiana, su producto también tiene una forma gaussiana, y luego la integral fea que escribimos anteriormente se toma fácilmente. La solución y su error serán los siguientes:

φ=(KTΣ1K+αΩ)1KTΣ1Tf

Σφ=(KTΣ1K+αΩ)1,

Dónde Σ- matriz de covarianza de la distribución gaussiana multidimensional, α- el valor más probable del parámetro α, que se determina a partir de la condición de máxima densidad de probabilidad a posteriori:

P(α|f)=CαRg(Ω)2|(KTΣ1K+αΩ)1|exp(12fTΣ1KT(KTΣ1K+αΩ)1KTΣ1Tf)


¿Y si no tengo errores gaussianos?


La segunda parte del artículo estará dedicada a esto, pero por ahora vamos a esbozar la esencia del problema.

φi=dφφiP(f|φ)dαP(α)αRg(Ω)2exp(α2(φ,Ωφ))dφP(f|φ)dαP(α)αRg(Ω)2exp(α2(φ,Ωφ))


El principal problema es que esta terrible integral, en primer lugar multidimensional, y en segundo lugar, en límites infinitos. Además, es un vector altamente multidimensionalφ puede tener dimensión fácilmente m=3050, y los métodos de cuadrícula para calcular integrales tienen una complejidad de tipo O(nm), por lo tanto, inaplicable en este caso. Al tomar integrales multidimensionales, la integración de Monte Carlo funciona bien.

Además, dado que nuestros límites son infinitos, debemos utilizar el método de muestreo importante, pero luego debemos seleccionar una función para el muestreo. Para automatizar todo, debe usar la Markov Chain Monte Carlo (MCMC) , que puede adaptar independientemente la función de muestreo a la integral calculada. Hablaremos sobre la aplicación de MCMC en el próximo artículo.

Parte práctica


La primera implementación del método de regularización estadística se escribió en los años 70 en Algol y se utilizó con éxito para los cálculos en física atmosférica. A pesar de que todavía tenemos las fuentes escritas a mano del algoritmo, decidimos agregar un poco de modernismo y hacer una implementación en Python, y luego en Julia.

Pitón


Instalación

Instalar a través de pip:

pip install statreg

o descargue el código fuente .

Ejemplos

Como ejemplo, considere cómo usar un módulo staregpara recuperar datos para una matriz y una ecuación integral.

Importamos los paquetes científicos necesarios.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.integrate import quad
%matplotlib inline

Determinamos la señal verdadera, que restauraremos.

a = 0
b = 5
#  
phi = lambda x: 4*norm.pdf(x-2, scale=0.4) + 2*norm.pdf(x-4, scale = 0.5)
x = np.linspace(a, b,100)
plt.plot(x, phi(x));


Defina el núcleo y la operación de convolución de funciones (Nota: np.convolutionespecíficamente para matrices):

kernel = lambda x,y : np.heaviside(x-y, 1) #  
convolution =  np.vectorize(lambda y: quad(lambda x: kernel(x,y)*phi(x), a,b)[0])

Generamos los datos medidos y los hacemos ruido usando la distribución normal:

y = np.linspace(a, b, 50)
ftrue = convolution(y)
sig = 0.05*ftrue +0.01 #  
f = norm.rvs(loc = ftrue, scale=sig)
plt.errorbar(y, f, yerr=sig);



Resolvemos la ecuación integral

Importamos el solucionador y la clase auxiliar para discretización:

from statreg.model import GaussErrorUnfolder
from statreg.basis import CubicSplines

Como base funcional para la discretización, usamos splines cúbicos y, como condición adicional, indicamos que la función toma valores cero en los bordes.

basis = CubicSplines(y, boundary='dirichlet')
model = GaussErrorUnfolder(basis, basis.omega(2))

Resolvemos la ecuación:

phi_reconstruct = model.solve(kernel, f, sig, y)

Estamos creando un cronograma de soluciones:

plt.plot(x,phi(x))
phir = phi_reconstruct(x)
phiEr = phi_reconstruct.error(x)
plt.plot(x, phir, 'g')
plt.fill_between(x, phir-phiEr, phir + phiEr, color='g', alpha=0.3);



Resolvemos la ecuación matricial

Importamos el solucionador y la clase auxiliar para discretización:

from statreg.model import GaussErrorMatrixUnfolder
from statreg.basis import CubicSplines

Para obtener matrices, utilizamos nuestra base funcional, pero está claro que las matrices se pueden obtener de cualquier manera.

cubicSplines = CubicSplines(y, boundary='dirichlet')
omega = cubicSplines.omega(2)
Kmn = cubicSplines.discretizeKernel(kernel,y)

Resolvemos la ecuación matricial:

model = GaussErrorMatrixUnfolder(omega)
result = model.solve(Kmn, f, sig)

Construye el gráfico:

phir = lambda x: sum([p*bf(x) for p, bf in zip(result.phi,cubicSplines.basisFun)])
plt.plot(x,phir(x))
plt.plot(x,phi(x));



Julia


Como mencionamos, un mayor desarrollo de la técnica requiere una integración avanzada de Monte Carlo. Podríamos usar algún módulo en Python (por ejemplo, trabajamos con PyMC3), sin embargo, estamos, entre otras cosas, participando en un proyecto conjunto con el Instituto Max Planck en Munich.

Este proyecto se llama Bayesian Analysis Toolkit . Su objetivo es crear un marco con herramientas para los métodos de análisis bayesianos, que incluyen principalmente herramientas para MCMC. Ahora el equipo está trabajando en la segunda versión del marco, que está escrita en Julia (la primera está escrita en C ++ incorrecto). Una de las tareas de nuestro grupo es demostrar las capacidades de este marco utilizando el ejemplo de regularización estadística, por lo que escribimos una implementación en Julia .

using PyCall
include("../src/gauss_error.jl")
include("../src/kernels.jl")

a = 0.
b = 6.

function phi(x::Float64)
    mu1 = 1.
    mu2 = 4.
    n1 = 4.
    n2 = 2.
    sig1 = 0.3
    sig2 = 0.5

    norm(n, mu, sig, x) = n / sqrt(2 * pi*sig^2) * exp(-(x - mu)^2 / (2 * sig^2))
    return norm(n1, mu1, sig1, x) + norm(n2, mu2, sig2, x)
end
x = collect(range(a, stop=b, length=300))

import PyPlot.plot

myplot = plot(x, phi.(x))
savefig("function.png", dpi=1000)


Esta vez usamos un núcleo diferente, no tomaremos un paso integrador, sino una convolución con un gaussiano, que en realidad conduce a un cierto "desenfoque" a nuestros datos:

function kernel(x::Float64, y::Float64)
    return getOpticsKernels("gaussian")(x, y)
end

convolution = y -> quadgk(x -> kernel(x,y) * phi(x), a, b, maxevals=10^7)[1]
y = collect(range(a, stop = b, length=50))
ftrue = convolution.(y)
sig = 0.05*abs.(ftrue) +[0.01 for i = 1:Base.length(ftrue)]
using Compat, Random, Distributions
noise = []
for sigma in sig
    n = rand(Normal(0., sigma), 1)[1]
    push!(noise, n)
end
f = ftrue + noise
plot(y, f)


Del mismo modo, tomamos la base de splines con extremos fijos:

basis = CubicSplineBasis(y, "dirichlet")
Kmn = discretize_kernel(basis, kernel, y)
model = GaussErrorMatrixUnfolder([omega(basis, 2)], "EmpiricalBayes", nothing, [1e-5], [1.], [0.5])
result = solve(model, Kmn, f, sig)
phivec = PhiVec(result, basis)

x = collect(range(a, stop=b, length=5000))
plot(x, phi.(x))

phi_reconstructed = phivec.phi_function.(x)
phi_reconstructed_errors = phivec.error_function.(x)

plot(x, phi_reconstructed)
fill_between(x, phi_reconstructed - phi_reconstructed_errors, phi_reconstructed + phi_reconstructed_errors, alpha=0.3)



Ejemplo del mundo real

Como ejemplo de análisis de datos reales, restauraremos el espectro de dispersión de electrones de una mezcla de hidrógeno-deuterio. En el experimento, se midió el espectro integral (es decir, el número de electrones está por encima de cierta energía), y necesitamos restaurar el espectro diferencial. Para estos datos, el espectro se reconstruyó inicialmente mediante el ajuste, por lo que tenemos una base para verificar la corrección de nuestro algoritmo.

Así es como se ve el espectro integrado inicial:

Y así, el resultado de la restauración: el

análisis con ajuste tiene tres desventajas principales:

  • , .
  • , .
  • .

La regularización estadística evita todos estos problemas y proporciona un resultado independiente del modelo con errores de medición. La solución obtenida por regularización está de acuerdo con la curva de ajuste. Tenga en cuenta los dos picos pequeños a 25 y 30 eV. Se sabe que se forma un pico a 25 eV durante la dispersión doble, y se restableció mediante un ajuste, ya que se especificó claramente en la función de ajuste. Un pico de 30 eV puede ser una anomalía estadística (los errores son bastante grandes en este punto) o, posiblemente, indica la presencia de dispersión disociativa adicional.

Conclusiones y anuncio de la próxima parte.


Le contamos sobre una técnica útil que se puede adaptar a muchas tareas de análisis de datos (incluido el aprendizaje automático) y obtener un "ajuste" honesto de la respuesta: la solución más racional a la ecuación ante la incertidumbre causada por errores de medición. Como un buen bono, obtenemos valores para el error de decisión. Aquellos que quieran participar en el desarrollo o aplicar el método de regularización estadística pueden contribuir en forma de código en Python, Julia o en otra cosa.

En la siguiente parte hablaremos sobre:

  • Usando MCMC
  • Descomposición de Cholesky
  • Como ejemplo práctico, consideramos el uso de la regularización para procesar una señal de un modelo de un detector orbital de protones y electrones.

Referencias



Publicado por Mikhail Zeleny , Investigador del Laboratorio de Métodos Experimentales de Física Nuclear en JetBrains Research .

All Articles