Aprendizaje anticausal en problemas inversos y su
aplicaci
´
on a tomograf
´
ıa optoac
´
ustica
Anticausal Learning for Inverse Problems and its Application on Optoacoustic Tomography
Matias Vera
∗†1
, Mart
´
ın G. Gonz
´
alez
∗†
, Leonardo Rey Vega
∗†
∗
Universidad de Buenos Aires, Facultad de Ingenier
´
ıa
Paseo Colon 850, C1063ACV, Buenos Aires, Argentina
†
Consejo Nacional de Investigaciones Cient
´
ıficas y T
´
ecnicas, (CONICET)
Godoy Cruz 2290, C1425FQB, Buenos Aires, Argentina
1
mvera@fi.uba.ar
Resumen— Los algoritmos de inteligencia artificial
habitualmente fallan cuando la distribuci
´
on de los datos
se desv
´
ıa de la utilizada durante el entrenamiento. Esta
vulnerabilidad puede ser corregida post-entrenamiento,
pero la misma puede requerir una etapa de ajuste
computacionalmente pesada y/o una gran necesidad de
nuevos datos. En este contexto, la teor
´
ıa de causalidad
suele ser un excelente paradigma para diferenciar los
mecanismos propensos a variaciones de los invariantes.
Esto permitir
´
ıa hacer un ajuste solamente sobre
el modelo variable, reduciendo la complejidad del
problema. Sin embargo, este paradigma est
´
a muy
poco estudiado en lo referido a problemas inversos,
principalmente porque estos problemas son por
definici
´
on anticausales. En este trabajo se analiza
el desempe
˜
no y limitaciones de algoritmos b
´
asicos
en problemas inversos que cumplan el requisito de
aprender de forma anticausal. En particular, se estudian
estos algoritmos en el contexto de reconstrucci
´
on de
im
´
agenes en tomograf
´
ıa optoac
´
ustica.
Palabras clave: problemas inversos; modelos guiados por
la f
´
ısica; teor
´
ıa de causalidad; tomograf
´
ıa optoac
´
ustica
Abstract— Artificial intelligence algorithms commonly
exhibit poor performance when deployed on data
whose distribution deviates from the one utilized during
the training phase. While this vulnerability can be
addressed post-training, doing so may necessitate a
computationally intensive fine-tuning process and/or
require a significant acquisition of new data. In this
context, causality theory presents an excellent paradigm
for distinguishing variation-prone mechanisms from
invariant ones. This distinction would permit fitting
the model exclusively to the variable components,
thereby reducing the complexity of the overall problem.
However, this paradigm remains under-explored in
relation to inverse problems, primarily because such
problems are, by their very definition, anticausal.
This work undertakes an analysis of the performance
and inherent limitations of fundamental algorithms
in inverse problems that satisfy the criteria for
anticausal learning. Specifically, these algorithms are
investigated within the context of image reconstruction
in optoacoustic tomography.
Keywords: inverse problems; physics-guided models;
causality theory; optoacoustic tomography
I. INTRODUCCI
´
ON
Los problemas inversos constituyen una clase particular
de tareas cuyo prop
´
osito es inferir causas desconocidas a
partir de efectos observados [1]. Se clasifican como tareas
anticausales porque las mismas buscan estimar causas a
partir de efectos, una direcci
´
on opuesta al proceso de ge-
neraci
´
on de las variables. Se presentan de manera natural
en numerosos contextos cient
´
ıficos y de ingenier
´
ıa, entre
ellos la reconstrucci
´
on de im
´
agenes, la exploraci
´
on geof
´
ısica
y el diagn
´
ostico m
´
edico. Ejemplos cl
´
asicos van desde la
reconstrucci
´
on de im
´
agenes de alta resoluci
´
on a partir de
mediciones degradadas en diagn
´
ostico por im
´
agenes [2] has-
ta la recuperaci
´
on de estructuras subterr
´
aneas mediante datos
s
´
ısmicos [3]. Estas tareas suelen ser mal condicionadas, es
decir, su soluci
´
on puede no ser
´
unica o extremadamente
sensible al ruido en las mediciones, lo que hace imprescin-
dible la incorporaci
´
on de conocimiento f
´
ısico previo para
caracterizar al problema. Se denominan entonces modelos
guiados por la f
´
ısica.
La tomograf
´
ıa optoac
´
ustica (TOA) es un m
´
etodo de
obtenci
´
on de im
´
agenes m
´
edicas mediante el uso del efecto
optoac
´
ustico. Un pulso de luz que incide en el tejido
biol
´
ogico blando se esparcir
´
a por el mismo y una parte ser
´
a
absorbida por mol
´
eculas presentes en la muestra biol
´
ogica,
conocidas como crom
´
oforos. La energ
´
ıa del crom
´
oforo
excitado se convierte luego en calor, que en el marco de
un proceso isoc
´
orico, termina generando un aumento de
presi
´
on [4]. Esto se detecta a trav
´
es de distintos arreglos de
sensores de ultrasonido, generando los llamados sinogramas:
representaciones gr
´
aficas de las se
˜
nales ac
´
usticas en funci
´
on
del tiempo medidas por cada detector. La gran cantidad
de configuraciones diferentes de medici
´
on en esta tarea,
as
´
ı como la presencia de incertidumbres o el conocimiento
parcial de los par
´
ametros, pueden dar lugar a algoritmos
de reconstrucci
´
on espec
´
ıficamente dise
˜
nados para una con-
Revista elektron, Vol. 9, No. 2, pp. 76-83 (2025)
ISSN 2525-0159
76
Recibido: 09/10/25; Aceptado: 09/11/25
https://doi.org/10.37537/rev.elektron.9.2.222.2025
Original Article
figuraci
´
on particular que podr
´
ıa no ser la que se utilizar
´
a
en una situaci
´
on pr
´
actica final, sufriendo un cambio en la
distribuci
´
on de los datos [5].
El Principio de Mecanismos Causales Independientes es
una hip
´
otesis heur
´
ıstica proveniente de la teor
´
ıa de cau-
salidad [6]. Este principio postula que el proceso causal
generativo de las variables de un sistema se compone de
m
´
odulos aut
´
onomos que no se informan ni influyen entre
s
´
ı. En un problema de TOA, esto se traduce en entrenar
por separado un modelo para la fuente (representados por
im
´
agenes) y otro para el proceso de medici
´
on de sinogramas
a partir de su fuente (problema directo). Una variaci
´
on
en la configuraci
´
on experimental solo afectar
´
ıa al segundo
modelo, el cu
´
al podr
´
ıa ser corregido con un proceso de adap-
taci
´
on de dominio [7] o aprendizaje por transferencia [8]. El
hecho de corregir solo una parte del modelo, podr
´
ıa aliviar
potencialmente el costo computacional y la gran necesidad
de datos del nuevo entorno. Este tipo de aprendizaje podr
´
ıa
incluso ayudar a definir algoritmos invariantes a cambios de
entorno [5], [9], [10].
El presente trabajo se centra en el estudio de algoritmos
b
´
asicos de aprendizaje anticausal en el marco de los pro-
blemas inversos y su aplicaci
´
on espec
´
ıfica a la tomograf
´
ıa
optoac
´
ustica. Con ello, se busca realizar una prueba de
concepto para evaluar la viabilidad y el potencial de la
convergencia de estos campos de estudio.
II. CONCEPTOS B
´
ASICOS
A. Problemas Inversos
A pesar de su generalidad, los problemas inversos siguen
un marco matem
´
atico bastante unificado. El objetivo es
recuperar una muestra desconocida y ∈ R
d
y
distribuida a
partir de p(y), suponiendo acceso a mediciones x ∈ R
d
x
y
asumiendo un modelo de la forma:
X = A(Y ) + V (1)
donde X, Y son variables aleatorias representativas de x
y y, y V ∼ N (0, Σ
v
) es un ruido aleatorio gaussiano
independiente de Y . En otras palabras, el modelo define
la relaci
´
on X|
Y =y
∼ N (A(y), Σ
v
). El predictor
´
optimo
en t
´
erminos de minimizar el error cuadr
´
atico medio es
E[Y |X = x], el cual se busca estimar:
E[Y |X = x] = arg min
f:R
d
x
→R
d
y
E[(f(X) − Y )
2
] (2)
Observaci
´
on 1: Muchas aplicaciones pr
´
acticas pue-
den aproximarse con un modelo lineal no invertible,
A(y) = A · y donde A es una matriz de dimensi
´
on
d
x
×d
y
. Esta matriz A suele estar fuertemente mal con-
dicionada, lo que vuelve escencial diferentes estrategias
de regularizaci
´
on para aproximar su inversi
´
on. En este
trabajo se tiene en cuenta esta hip
´
otesis, entonces se
obtiene el siguiente modelo directo de (1):
x = Ay + v (3)
B. Fuera de Distribuci
´
on
En aplicaciones de sensado es razonable suponer que la
distribuci
´
on de los datos puede cambiar porque las condi-
ciones de adquisici
´
on y los factores experimentales rara vez
son id
´
enticos entre mediciones. Por ejemplo, en im
´
agenes
de TOA, la posici
´
on del sensor, la velocidad del sonido o
incluso peque
˜
nas variaciones en la muestra pueden alterar
la relaci
´
on entre la fuente imagen original) y la observaci
´
on
x (sinograma).
Sea p
e
(x, y) la distribuci
´
on de los datos, denotamos e ∈ E
a las posibles variaciones que puede sufrir. La cuesti
´
on
clave a abordar es la representaci
´
on de los entornos en el
contexto de un problema inverso. Un ejemplo paradigm
´
atico
de problema inverso en ingenier
´
ıa involucra un escenario
en el cual y representa una variable f
´
ısica a sensar, y x
corresponde a su medici
´
on indirecta. Si se supone que las
posibles variaciones se deben a cambios en las condiciones
experimentales donde se lleva a cabo el sensado, es razona-
ble suponer que la distribuci
´
on p(y) permanece fija mientras
que el componente dependiente del entorno queda capturado
por p
e
(x|y).
En la teor
´
ıa de la causalidad, este tipo de problema se
conoce como aprendizaje anticausal, dado que la meta es
predecir la causa a partir del efecto. En estos problemas,
el enfoque recomendado a seguir es el que se describe a
continuaci
´
on [6]:
Observaci
´
on 2: p
e
(x|y) representa el mecanismo
causal que genera X a partir de Y , y es independiente
de la distribuci
´
on de la causa, p(y). Por otro lado,
p
e
(y|x) es sensible al cambio en la distribuci
´
on de
p(y). Por lo tanto, en t
´
erminos generales, al estimar
p
e
(y|x) conviene modelar p
e
(x|y) y p(y) por separado
y luego construir p
e
(y|x) usando la regla de Bayes.
En este contexto, se estudian diversas t
´
ecnicas para
estimar E[Y |X = x], con el objetivo de desacoplar el
modelado de p(x|y) de p(y), aprovechando el conocimiento
de la distribuci
´
on X|
Y =y
∼ N (Ay, Σ
v
). La estimaci
´
on se
efectuar
´
a en dos etapas:
1. Modelo directo: estimar los par
´
ametros (A, Σ
v
) a
partir de un conjunto de datos aleatorio {(X
i
, Y
i
)}
n
i=1
.
En este trabajo, para el modelo directo se supondr
´
a:
1) una matriz de covarianza esf
´
erica para el ruido
Σ
v
= σ
2
v
· I, donde σ
2
v
es la estimaci
´
on emp
´
ırica de
la varianza; y 2) la matriz A ser
´
a determinada con la
metodolog
´
ıa est
´
andar de TOA (la cu
´
al ser
´
a explicada
a continuaci
´
on).
2. Componente invariante: desarrollar un procedimien-
to capaz de realizar este modelado desacoplado de
p(x|y) (conocido en este paso) y p(y) (desconocido),
sin comprometer significativamente el rendimiento en
comparaci
´
on con los m
´
etodos cl
´
asicos (dise
˜
nados para
un
´
unico entorno). En esta etapa se utiliza un conjunto
de datos simple de objetivos {Y
i
}
n
i=1
. Nos enfoca-
remos en algoritmos no neuronales, con el fin de
compararlos con m
´
etodos cl
´
asicos como la proyecci
´
on
lineal inversa (LBP, por sus siglas en ingl
´
es).
Revista elektron, Vol. 9, No. 2, pp. 76-83 (2025)
ISSN 2525-0159
77
http://elektron.fi.uba.ar
C. Tomograf
´
ıa Optoac
´
ustica
Es bien sabido que, tras la excitaci
´
on de una muestra
biol
´
ogica por un pulso electromagn
´
etico δ(t), la presi
´
on
ac
´
ustica p(r, t) en la posici
´
on r ∈ R
3
y tiempo t satisface
la ecuaci
´
on diferencial [11]:
∂
2
∂t
2
− v
2
s
∇
2
p(r, t) = 0 (4)
con las condiciones iniciales:
p(r, 0) = p
0
(r) ,
∂p
∂t
(r, 0) = 0 (5)
donde p
0
(r) es la presi
´
on optoac
´
ustica inicial y v
s
representa
la velocidad del sonido en el medio, el cual se supone
homog
´
eneo y sin absorci
´
on ac
´
ustica. Bajo la hip
´
otesis
usual de confinamiento t
´
ermico y ac
´
ustico [12], es decir,
cuando la duraci
´
on del pulso l
´
aser es lo suficientemente
corta como para que se pueda despreciar la conducci
´
on de
calor y la propagaci
´
on ac
´
ustica hacia regiones vecinas de
la zona iluminada, la presi
´
on inducida inicialmente p
0
(r) es
proporcional a la densidad de energ
´
ıa
´
optica total absorbida.
Usando la formalizaci
´
on de la funci
´
on de Green, la presi
´
on
recibida por un detector puntual ideal en la posici
´
on r
d
puede escribirse como:
p
d
(r
d
, t) =
1
4π v
2
s
∂
∂t
ZZZ
V
p
0
(r)
δ (t − |r
d
− r|/v
s
)
|r
d
− r|
dr
(6)
El objetivo del problema inverso en TOA es reconstruir
p
0
(r) a partir del sinograma p
d
(r
d
, t) medido en varias
posiciones r
d
, que t
´
ıpicamente se encuentran sobre una
superficie S que contiene el volumen de inter
´
es [13].
Varios enfoques, como los algoritmos de retroproyecci
´
on
[14], [15], son de los m
´
as populares y utilizados en el
problema de reconstrucci
´
on de im
´
agenes en TOA. Dichos
m
´
etodos proporcionan f
´
ormulas de reconstrucci
´
on en forma
cerrada en t
´
erminos de las se
˜
nales detectadas sobre la
superficie de detecci
´
on. Sin embargo, estos m
´
etodos suponen
que los detectores son puntuales, sin limitaciones de ancho
de banda y con respuesta angular isotr
´
opica [16]. En la
pr
´
actica, los transductores tienen tama
˜
no finito, ancho de
banda limitado y su respuesta espacial no es constante [17].
Adem
´
as, las se
˜
nales detectadas son ruidosas. Estas desvia-
ciones del escenario ideal supuesto por las f
´
ormulas exactas
de reconstrucci
´
on pueden generar artefactos e im
´
agenes
distorsionadas.
Un enfoque diferente pero relacionado al problema de
reconstrucci
´
on est
´
a dado por los algoritmos basados en
matrices [18]. En esta t
´
ecnica, la soluci
´
on directa de (6)
se discretiza. Como resultado, se obtiene la ecuaci
´
on ma-
tricial (3) que se usa para resolver el problema inverso,
donde x es un vector columna que representa las presiones
medidas en un conjunto de posiciones de detectores r
d
l
(l = 1 . . . N
d
) y en instantes de tiempo t
k
(k = 1 . . . N
t
);
y es un vector columna que representa los valores de la
presi
´
on ac
´
ustica inicial, y que t
´
ıpicamente se denomina la
imagen de referencia; y A es la matriz del modelo. El
j-
´
esimo elemento (j = 1 . . . N) en y contiene el valor
promedio de la presi
´
on inicial dentro de un elemento de
volumen de tama
˜
no ∆V en la posici
´
on r
j
. Una de las
ventajas de este enfoque es que cualquier efecto lineal en
el sistema puede ser considerado f
´
acilmente (por ejemplo,
la respuesta espacial y temporal de los sensores) [19]. Una
vez establecida la formulaci
´
on discreta, el problema inverso
se reduce al problema algebraico de invertir (3). La matriz
A puede escribirse como la multiplicaci
´
on de dos matrices
A
oa
A
s
, donde A
s
representa la funci
´
on de respuesta del
sistema de imagen para un sensor puntual ideal y A
oa
es
la forma matricial de un operador de derivada temporal. La
matriz A
s
se define como [20]:
A
s
lkj
=
1
4πv
2
s
∆V
∆t
2
d(t
k
, r
j
, r
dl
)
|r
d
l
− r
j
|
(7)
d(t
k
, r
j
, r
dl
) =
(
1 si |t
k
−
|r
d
l
−r
j
|
v
s
| < ∆t/2
0 en otro caso
(8)
donde ∆t es el paso temporal en el que se muestrean las
se
˜
nales p
d
(r
d
, t). No es dif
´
ıcil ver que (7) constituye una
discretizaci
´
on del integrando en (6), mientras que (8) indica
el instante de tiempo en que el efecto de la presi
´
on inicial en
la posici
´
on r
j
es capturado por el sensor r
d
l
. En el caso de
un detector de tama
˜
no finito, la respuesta impulsiva espacial
del sensor se tiene en cuenta dividiendo el
´
area del sensor en
elementos superficiales (tratados como detectores puntuales)
que luego se suman [20], [21].
Observaci
´
on 3: Es razonable suponer que, en una
aplicaci
´
on real, la matriz del modelo no se conoce
completamente. En este trabajo, consideraremos que
la diferencia entre ambas reside en la incertidumbre de
la posici
´
on del sensor y de la velocidad del sonido al
calcularla.
III. MODELOS ANTICAUSALES
Nos enfocamos en modelos que definen una funci
´
on de
aproximaci
´
on ˆy(x) basada en los par
´
ametros (A, σ
2
v
) y un
conjunto de datos objetivo {Y
i
}
n
i=1
. Estos modelos buscan
abordar los desaf
´
ıos de los problemas inversos aprovechan-
do representaciones anticausales desacopladas y estrategias
computacionales eficientes.
A. Proyecci
´
on Lineal Inversa (LBP)
El problema de inversi
´
on puede formularse usando un cri-
terio cuadr
´
atico combinado con un t
´
ermino de regularizaci
´
on
de Tikhonov:
ˆy
TIK
(x) = m´ın
y∈R
d
y
∥Ay − x∥
2
+ λ∥y∥
2
(9)
donde λ ≥ 0 es un par
´
ametro de regularizaci
´
on que
mejora la estabilidad del problema inverso, el cual suele ser
mal condicionado. Adem
´
as, este t
´
ermino de regularizaci
´
on
mitiga los efectos del ruido en las se
˜
nales medidas. La solu-
ci
´
on al problema regularizado de Tikhonov puede derivarse
anal
´
ıticamente como [22]:
ˆy
TIK
(x) = (A
T
A + λI)
−1
A
T
x (10)
Sin embargo, calcular esta inversa puede ser compu-
tacionalmente costoso, particularmente en problemas a gran
escala. Para abordar esto, una simplificaci
´
on com
´
un es
Revista elektron, Vol. 9, No. 2, pp. 76-83 (2025)
ISSN 2525-0159
78
http://elektron.fi.uba.ar
considerar valores grandes de λ. En el r
´
egimen asint
´
oti-
co cuando λ → ∞, la soluci
´
on de Tikhonov se vuelve
proporcional al operador adjunto ˆy
TIK
(x) ∝ A
T
x, y esta
aproximaci
´
on conduce a una soluci
´
on computacionalmente
eficiente conocida como proyecci
´
on lineal inversa (LBP)
[23]:
ˆy
LBP
(x) = A
T
x (11)
La soluci
´
on LBP proporciona un m
´
etodo sencillo y
computacionalmente econ
´
omico para obtener una recons-
trucci
´
on inicial de la imagen. Sin embargo, es importante
notar que las reconstrucciones basadas en LBP suelen exhi-
bir limitaciones y artefactos, particularmente en escenarios
de visi
´
on limitada [24]. A pesar de estas desventajas, LBP
sigue siendo ampliamente utilizada debido a su simplicidad
y facilidad de implementaci
´
on. Cabe destacar que esta apro-
ximaci
´
on depende
´
unicamente de la matriz A, independiente
de σ
2
v
y del conjunto de datos.
B. Estimaci
´
on Monte Carlo
El c
´
alculo de E[Y |X = x] mediante la descomposici
´
on
causal p(x, y) = p(y)p(x|y) para este modelo puede escri-
birse como:
E[Y |X = x] =
Z
R
d
y
y
p(x|y)p(y)
p(x)
dy (12)
=
R
R
d
y
ye
−
1
2
(x−Ay)
T
Σ
−1
v
(x−Ay)
p(y)dy
R
R
d
y
e
−
1
2
(x−Ay)
T
Σ
−1
v
(x−Ay)
p(y)dy
(13)
Una estimaci
´
on Monte Carlo puede realizarse de manera
directa como:
ˆy
MC
(x) =
P
n
i=1
y
i
e
−
1
2
(x−Ay
i
)
T
Σ
−1
v
(x−Ay
i
)
P
n
i=1
e
−
1
2
(x−Ay
i
)
T
Σ
−1
v
(x−Ay
i
)
C. Modelo de Mezcla Gaussiana (GMM)
El enfoque propuesto consiste en aprender p(y) median-
te un Modelo de Mezcla Gaussiana (GMM) entrenado a
trav
´
es del algoritmo de Expectativa-Maximizaci
´
on (EM) con
m gaussianas diagonales. Para integrar este modelo con
X|Y = y, se emplea un marco de an
´
alisis de factores.
El GMM supone la existencia de una variable mezcladora
categ
´
orica K ∼ Cat(ω
1
, · · · , ω
m
) cuyo v
´
ınculo con el resto
de las variables est
´
a definido por la cadena de Markov
K → Y → X (donde las flechas definen la relaci
´
on de
causalidad supuesta). Es decir, el modelo supone por un
lado la distribuci
´
on Y |K = k ∼ N (µ
k
, Λ
k
) (con Λ
k
una
matriz diagonal) y por el otro X|Y = y ∼ N (Ay, σ
2
v
· I).
El algoritmo EM entrena {ω
k
, µ
k
, Λ
k
}
m
k=1
a partir de
un conjunto de datos {Y
i
}
n
i=1
. Tras el entrenamiento, la
inferencia puede realizarse con el siguiente lema.
Lema 1: En el modelo descrito anteriormente, la
esperanza condicional puede calcularse como:
E[Y |X = x] =
m
X
k=1
P(K = k|X = x)E[Y |X = x, K = k]
(14)
donde
E[Y |X = x, K = k] (15)
= µ
k
+ Λ
k
A
T
σ
2
v
I + AΛ
k
A
T
−1
(x − Aµ
k
)
y
P(K = k|X = x) =
ω
k
· N
x
Aµ
k
, σ
2
v
I + AΛ
k
A
T
P
m
l=1
ω
l
· N
x
(Aµ
l
, σ
2
v
I + AΛ
l
A
T
)
(16)
con N
x
(µ, Σ) la densidad de probabilidad asociada a
una distribuci
´
on normal de media µ y covarianza Σ
evaluada en x.
La demostraci
´
on puede verse en el Ap
´
endice A. Dentro
de este marco, definimos el estimador ˆy
GMM
(x) = E[Y |X =
x], calculado usando el Lema 1.
D. F
´
ormula de Tweedie
El operador de esperanza condicional es una herramienta
fundamental en numerosos campos que dependen del an
´
ali-
sis estad
´
ıstico. Existen diversas identidades derivadas que
establecen relaciones entre la esperanza condicional y otras
cantidades estad
´
ısticas, como la varianza condicional. Entre
estas, una de las identidades m
´
as relevantes es la f
´
ormula
de Tweedie [25], que proporciona un m
´
etodo para calcular
la esperanza condicional a trav
´
es de la medida marginal.
En este contexto, la f
´
ormula de Tweedie se formaliza en el
siguiente lema.
Lema 2: La esperanza condicional E[Y |X = x]
satisface la siguiente identidad:
A · E[Y |X = x] = x + σ
2
v
· ψ(x) (17)
donde ψ(x) = ∇
x
log p(x) denota la funci
´
on score,
con ∇
x
representando el gradiente respecto de x. La
funci
´
on score tambi
´
en puede expresarse como:
ψ(x) = ∇
x
log E
e
−
1
2σ
2
v
∥x−AY ∥
2
(18)
La demostraci
´
on se encuentra en el Ap
´
endice B. En este
contexto, definimos la aproximaci
´
on de Tweedie como:
ˆy
TW
(x) = A
T
x + σ
2
v
A
T
· ψ(x) (19)
Dentro de este marco, la funci
´
on de Tweedie puede verse
como una generalizaci
´
on de LBP (11), donde la desviaci
´
on
respecto de la formulaci
´
on est
´
andar de LBP est
´
a gobernada
por el par
´
ametro σ
2
v
, estimado en el modelo directo. La apro-
ximaci
´
on de la funci
´
on score, sin embargo, sigue siendo un
problema abierto. Nosotros utilizaremos una aproximaci
´
on
de Taylor de primer orden de la siguiente manera:
e
−
1
2σ
2
v
∥x−Ay∥
2
≈ 1 −
1
2σ
2
v
(x
T
x − 2y
T
A
T
x + y
T
A
T
Ay)
= 1 −
1
2σ
2
v
(x
T
x − 2y
T
A
T
x + Tr(A
T
Ayy
T
)) (20)
donde Tr(·) es la traza de la matriz. La aproximaci
´
on de la
Revista elektron, Vol. 9, No. 2, pp. 76-83 (2025)
ISSN 2525-0159
79
http://elektron.fi.uba.ar
funci
´
on score se da entonces por:
ψ(x) =
∇
x
E
e
−
1
2σ
2
v
∥x−AY ∥
2
E
h
e
−
1
2σ
2
v
∥x−AY ∥
2
i
≈
−
1
σ
2
v
(x − AE[Y ])
1 −
1
2σ
2
v
(x
T
x − 2E[Y ]
T
A
T
x + Tr(A
T
AE[Y Y
T
])
(21)
Definimos la funci
´
on de aproximaci
´
on de Taylor como:
ˆy
TAY
(x) =
A
T
x −
(A
T
x − A
T
AE[Y ])
1 −
1
2σ
2
v
(x
T
x − 2E[Y ]
T
A
T
x + Tr(A
T
AE[Y Y
T
])
(22)
Cuando σ
2
v
→ 0, esta aproximaci
´
on converge a la soluci
´
on
LBP ˆy
TAY
(x) = ˆy
LBP
(x) = A
T
x, pero, por el contrario,
cuando σ
2
v
→ ∞, la aproximaci
´
on se vuelve constante
ˆy
TAY
(x) = A
T
AE[Y ]. En este contexto, σ
2
v
sirve como
medida de la calidad de la matriz A: si el ruido es peque
˜
no,
la reconstrucci
´
on se basa principalmente en A; si el ruido
es significativo, la reconstrucci
´
on tambi
´
en incorpora infor-
maci
´
on del conjunto de datos a trav
´
es de E[Y ] y E[Y Y
T
]
(estimados de forma emp
´
ırica).
IV. RESULTADOS EXPERIMENTALES
Los siguientes ejemplos se realizaron en una aplicaci
´
on
de TOA, donde x es el sinograma e y la imagen original.
Utilizamos un conjunto de entrenamiento de 5652 im
´
agenes
{y
i
}
5652
i=1
. Dicho conjunto est
´
a conformado por datos sint
´
eti-
cos generados a partir de im
´
agenes de vasculatura retiniana
obtenidas de bases de datos p
´
ublicas [26]–[30]. Tanto la
matriz A como la naturaleza del ruido aditivo se conocer
´
an
parcialmente. Las incertidumbres en la posici
´
on del sensor
y la velocidad del sonido en la matriz A son del 0,5 %.
Las im
´
agenes son creadas con ruido blanco gaussiano de
varianza aleatoria, elegida para que el valor de la SNR del
sinograma var
´
ıe uniformemente entre 20 dB y 80 dB (σ
2
v
ser
´
a estimado emp
´
ıricamente y se considerar
´
a fijo).
Por
´
ultimo, se aplica aumento de datos para incrementar
el conjunto de entrenamiento a 50000 muestras utilizando
datos sint
´
eticos generados por el modelo de difusi
´
on pre-
sentado en [31].
A. Estimaci
´
on Monte-Carlo
El primer experimento utiliza la estimaci
´
on Monte-Carlo
ˆy
MC
(x). Los resultados fueron interesantes: aproximadamen-
te el 70 % de las estimaciones resultaron ser solo ruido,
mientras que el 30 % restante tuvo un desempe
˜
no muy
bueno. Los resultados exitosos se muestran en la Fig. 1. Al
aplicar aumento de datos en este experimento, los resultados
fueron desalentadores: ninguna de las reconstrucciones fue
exitosa, observ
´
andose ruido en todos los casos. Se concluye
entonces que este m
´
etodo es muy sensible a la falta de
muestras en un entorno del punto de inferencia. El espacio
imagen posee una complejidad inherente tal, que la cantidad
de las muestras con las que se cuenta no alcanza para hacer
converger a la ley de los grandes n
´
umeros.
B. Modelo de Mezcla de Gaussianas
El segundo experimento utiliza el GMM previamente
descrito ˆy
GMM
(x) con m = 100. Una verificaci
´
on inicial
de la calidad de aprendizaje del GMM es usarlo como
m
´
etodo generativo. Sin aumento de datos, las im
´
agenes
generadas fueron muy ruidosas y presentaban pocas venas.
Sin embargo, con aumento de datos, el desempe
˜
no fue
satisfactorio. Los ejemplos se presentan en la Fig. 2. Aunque
las im
´
agenes generadas parecen razonables, el desempe
˜
no
del GMM en reconstrucci
´
on de im
´
agenes no fue satisfac-
toria, produciendo salidas ruidosas como puede verse en la
segunda fila de la Fig. 3.
Un an
´
alisis de los par
´
ametros optimizados revela que
las gaussianas exhiben covarianzas relativamente estrechas.
Este resultado implica que la regi
´
on del espacio latente
correctamente caracterizada se concentra principalmente en
torno a las m medias de las gaussianas. Este fen
´
omeno
podr
´
ıa atribuirse a que el aumento de datos sesg
´
o el proceso
de aprendizaje hacia las muestras reales disponibles. Es
plausible que la estrategia de aumento de datos empleada no
sea la m
´
as adecuada para la estimaci
´
on anticausal propuesta,
considerando, adem
´
as, su efecto perjudicial previamente
observado en la evaluaci
´
on mediante Monte Carlo. No
obstante, dado que el algoritmo no lograba converger a los
patrones esperados sin la aplicaci
´
on del aumento de datos,
se infiere una doble limitaci
´
on: por un lado, el n
´
umero
de muestras disponibles es insuficiente para el n
´
umero
de gaussianas utilizadas (justificando el uso del aumento
de datos); por otro lado, la complejidad intr
´
ınseca de las
im
´
agenes impide una reducci
´
on del tama
˜
no del conjunto
de datos. En esencia, nos enfrentamos a un compromiso
metodol
´
ogico o una limitaci
´
on fundamental del conjunto de
datos.
C. F
´
ormula de Tweedie
El tercer experimento utiliza la aproximaci
´
on de Taylor
de la f
´
ormula de Tweedie ˆy
TAY
(x). La Fig. 3 compara
el desempe
˜
no de este m
´
etodo con la aproximaci
´
on LBP,
mostrando resultados similares. En la Fig. 4, se puede
observar el desempe
˜
no de ˆy
TAY
(x) en funci
´
on de σ
2
v
para
diferentes m
´
etricas populares:
´
Indice de Similitud Estructu-
ral (SSIM), Correlaci
´
on de Pearson (PC), Error Cuadr
´
atico
Medio (RMSE) y Relaci
´
on Se
˜
nal-Ruido de Pico (PSNR) [5].
Mientras que SSIM, RMSE y PSNR son
´
optimos cuando
σ
2
v
= 0 (aproximaci
´
on LBP), en t
´
erminos de PC existe una
destacable ventaja al usar la aproximaci
´
on de Taylor.
V. CONCLUSIONES Y TRABAJOS FUTUROS
En el presente trabajo se propuso y se evalu
´
o una me-
todolog
´
ıa para la resoluci
´
on de problemas inversos desde
una perspectiva anticausal, la cual facilita la adaptaci
´
on a
posteriores cambios en la distribuci
´
on de los datos (cambios
de entorno). Se exploraron m
´
etodos b
´
asicos, incluyendo
LBP, Monte Carlo, GMM y Tweedie. Si bien LBP demostr
´
o
ser el m
´
etodo m
´
as consistentemente robusto, fue superado en
rendimiento en escenarios espec
´
ıficos. El m
´
etodo de Monte
Carlo mostr
´
o la capacidad de alcanzar una alta precisi
´
on en
regiones puntuales del espacio de soluci
´
on, y el estimador
de Tweedie evidenci
´
o mejoras significativas con respecto a
m
´
etrica PC.
Revista elektron, Vol. 9, No. 2, pp. 76-83 (2025)
ISSN 2525-0159
80
http://elektron.fi.uba.ar
Figura 1: Ejemplos de las im
´
agenes reconstruidas con Monte-Carlo para los casos donde el m
´
etodo funciona. La primera
fila muestra la imagen original y la segunda la reconstrucci
´
on Monte-Carlo.
Figura 2: Desempe
˜
no del GMM con aumento de datos como modelo generativo.
Figura 3: Ejemplos de im
´
agenes reconstruidas. La primera fila muestra la imagen original, la segunda la reconstrucci
´
on
LBP, la tercera GMM (algoritmo EM) entrenada con aumento de datos, y la cuarta utilizando aproximaci
´
on de Taylor de
la f
´
ormula de Tweedie.
La descomposici
´
on anticausal p(x|y) − p(y), presenta
muchas potencialidades una vez aprendidas, pero evidencia
un aumento en la dificultad de dicho aprendizaje. Mientras
que todos los m
´
etodos utilizaron el mismo modelado de
p(x|y) (caracterizados por A y σ
2
v
), los de mejor desempe
˜
no
fueron LBP (el cu
´
al evita modelar p(y)) y el estimador
Tweedie (solamente modela los momentos E[Y ] y E[Y Y
T
]),
dando a entender que la dificultad est
´
a en el aprendizaje del
espacio de las im
´
agenes. Esto se corrobora al analizar el
pobre desempe
˜
no del aumento de datos utilizado, el cual
ya hab
´
ıa demostrado buen desempe
˜
no en configuraciones
est
´
andar [31]. Es evidente entonces que la complejidad de
aprendizaje en el espacio de las im
´
agenes es superior al del
espacio de sinogramas.
El an
´
alisis aqu
´
ı presentado establece una base inicial en
esta
´
area de estudio, dejando un amplio espectro para la
investigaci
´
on futura. Una l
´
ınea de exploraci
´
on pendiente
consiste en optimizar la aproximaci
´
on de Tweedie para ga-
rantizar una utilizaci
´
on m
´
as eficiente de los datos, superan-
do la informaci
´
on provista
´
unicamente por los momentos
estad
´
ısticos. Adicionalmente, ser
´
ıa de gran relevancia incor-
porar una m
´
etrica de desempe
˜
no que est
´
e directamente aco-
Revista elektron, Vol. 9, No. 2, pp. 76-83 (2025)
ISSN 2525-0159
81
http://elektron.fi.uba.ar
Figura 4: Figuras de m
´
erito (SSIM, PC, RMSE y PSNR) en funci
´
on de σ
2
v
usando la aproximaci
´
on de Taylor.
plada con la adaptaci
´
on efectiva al entorno (trascendiendo la
mera potencialidad te
´
orica) para as
´
ı cuantificar su beneficio
real. Finalmente, resulta crucial investigar la adaptaci
´
on de
este enfoque anticausal a arquitecturas avanzadas, como las
redes neuronales o los modelos difusivos. Este esfuerzo
implicar
´
ıa no solo la modelizaci
´
on del problema directo,
sino tambi
´
en la incorporaci
´
on expl
´
ıcita de la identidad de
Bayes, como se demostr
´
o en trabajos previos (por ejemplo,
(13)). Sin embargo, la definici
´
on de la estrategia de imple-
mentaci
´
on
´
optima requiere exploraci
´
on adicional.
AGRADECIMIENTOS
Este trabajo fue financiado por la Universidad de Bue-
nos Aires (UBACYT 20020190100032BA), CONICET (PIP
11220200101826CO) y la Agencia I+D+i (PICT 2020-
01336).
REFERENCIAS
[1] G. Daras, H. Chung, C.-H. Lai, Y. Mitsufuji, J. C. Ye, P. Milanfar,
A. Dimakis, and M. Delbracio, “A survey on diffusion models for
inverse problems,” CoRR, 2024, arXiv e-prints: 2410.00083.
[2] Y. Song, L. Shen, L. Xing, and S. Ermon, “Solving inverse problems
in medical imaging with score-based generative models,” in Interna-
tional Conference on Learning Representations (ICLR), 2022.
[3] P. Lailly, “The seismic inverse problem as a sequence of before
stack migrations,” in Conference on Inverse Scattering, Theory and
application, 1983, pp. 206–220.
[4] M. Xu and L. Wang, “Photoacoustic imaging in biomedicine,” Rev.
Sci. Instrum., vol. 77, p. 041101, 2006.
[5] M. Vera, M. G. Gonz
´
alez, and L. R. Vega, “Invariant representations
in deep learning for optoacoustic imaging,” Review of Scientific
Instruments, vol. 94, no. 5, p. 054904, 05 2023. [Online]. Available:
https://doi.org/10.1063/5.0139286
[6] B. Sch
¨
olkopf, D. Janzing, J. Peters, E. Sgouritsa, K. Zhang, and
J. Mooij, “On causal and anticausal learning,” in International Con-
ference on Machine Learning (ICML), ser. ICML’12. Madison, WI,
USA: Omnipress, 2012, p. 459–466.
[7] H. Daum
´
e III, “Frustratingly easy domain adaptation,” in Annual
Meeting of the Association of Computational Linguistics (ACL).
Prague, Czech Republic: Association for Computational Linguistics,
2007, pp. 256–263.
[8] M. Rojas-Carulla, B. Sch
¨
olkopf, R. Turner, and J. Peters, “Invariant
models for causal transfer learning,” Journal of Machine Learning
Research (JMLR), vol. 19, no. 36, pp. 1–34, 2018.
[9] M. Arjovsky, L. Bottou, I. Gulrajani, and D. Lopez-Paz, “Invariant
risk minimization,” CoRR, Jul. 2019, arXiv e-prints: 1907.02893.
[10] B. Aubin, A. Slowik, M. Arjovsky, L. Bottou, and D. Lopez-Paz,
“Linear unit-tests for invariance discovery,” CoRR, Feb. 2021, arXiv
e-prints: 2102.10867.
[11] L. Wang and H. Wu, Biomedical Optics: Principles and Imaging,
1st ed. Hoboken, N.J: Wiley-Interscience, May 2007.
[12] R. Kruger, P. Liu, Y. Fang, and R. Appledorn, “Photoacoustic ultra-
sound (paus)—reconstruction tomography,” Medical Physics, vol. 22,
no. 10, p. 1605–1609, 1995.
[13] C. Lutzweiler and D. Razansky, “Optoacoustic imaging and to-
mography: reconstruction approaches and outstanding challenges in
image performance and quantification,” Sensors, vol. 13, pp. 7345–
7384, 2013.
[14] M. Xu and L. Wang, “Universal back-projection algorithm for pho-
Revista elektron, Vol. 9, No. 2, pp. 76-83 (2025)
ISSN 2525-0159
82
http://elektron.fi.uba.ar
toacoustic computed tomography,” Phys. Rev. E, vol. 71, p. 016706,
2005.
[15] P. Burgholzer, J. Bauer-Marschallinger, H. Gr
¨
un, M. Haltmeier, and
G. Paltauf, “Temporal back-projection algorithms for photoacous-
tic tomography with integrating line detectors,” Inverse Problems,
vol. 23, no. 6, p. S65–S80, Nov 2007.
[16] A. Rosenthal, V. Ntziachristos, and D. Razansky, “Acoustic inversion
in optoacoustic tomography: A review,” Current Medical Imaging
Reviews, vol. 9, pp. 318–336, 2013.
[17] C. Tian, M. Pei, K. Shen, S. Liu, Z. Hu, and T. Feng, “Impact
of system factors on the performance of photoacoustic tomography
scanners,” Phys. Rev. Applied, vol. 13, p. 014001, 2020.
[18] A. Rosenthal, D. Razansky, and V. Ntziachristos, “Fast semi-analytical
model-based acoustic inversion for quantitative optoacoustic tomo-
graphy,” IEEE Transactions on Medical Imaging, vol. 29, no. 6, p.
1275–1285, Jun 2010.
[19] L. Hirsch, M. G. Gonz
´
alez, and L. Rey Vega, “On the robustness of
model-based algorithms for photoacoustic tomography: Comparison
between time and frequency domains,” Review of Scientific Instru-
ments, vol. 92, no. 11, p. 114901, Nov 2021.
[20] G. Paltauf, P. R. Torke, and R. Nuster, “Modeling photoacoustic ima-
ging with a scanning focused detector using monte carlo simulation
of energy deposition,” Journal of biomedical optics, vol. 23, no. 12,
2018.
[21] A. Rosenthal, V. Ntziachristos, and D. Razansky, “Model-based op-
toacoustic inversion with arbitrary-shape detectors,” Medical Physics,
vol. 38, no. 7, pp. 4285–4295, 2011.
[22] A. D. Cezaro and F. T. D. Cezaro, “Regularization approaches for
quantitative photoacoustic tomography using the radiative transfer
equation,” Journal of Computational and Applied Mathematics, vol.
292, pp. 1–14, 2015.
[23] S. R. Arridge, M. M. Betcke, B. T. Cox, F. Lucka, and B. E. Treeby,
“On the adjoint operator in photoacoustic tomography,” Inverse Pro-
blems, vol. 32, no. 11, p. 115012, Oct 2016.
[24] M. G. Gonz
´
alez, M. Vera, and L. J. R. Vega, “Combining band-
frequency separation and deep neural networks for optoacoustic
imaging,” Optics and Lasers in Engineering, vol. 163, p. 107471,
2023.
[25] A. Hyv
¨
arinen, “Estimation of non-normalized statistical models by
score matching,” Journal of Machine Learning Research, vol. 6,
no. 24, pp. 695–709, 2005.
[26] DRIVE, “DRIVE: Digital retinal images for vessel extraction,” 2020.
[Online]. Available: https://drive.grand-challenge.org/
[27] ARIA, “Automated retinal image analysis,” 2006. [Online]. Available:
http://www.damianjjfarnell.com/
[28] RITE, “Retinal images vessel tree extraction,” 2013. [Online].
Available: https://medicine.uiowa.edu/eye/rite-dataset
[29] STARE, “Structured analysis of the retina,” 2000. [Online]. Available:
https://cecas.clemson.edu/
∼
ahoover/stare/
[30] A. Hatamizadeh, H. Hosseini, N. Patel, J. Choi, C. Pole, C. Hoeferlin,
S. Schwartz, and D. Terzopoulos, “RAVIR: A dataset and methodo-
logy for the semantic segmentation and quantitative analysis of retinal
arteries and veins in infrared reflectance imaging,” IEEE Journal of
Biomedical and Health Informatics, 2022.
[31] M. Gonz
´
alez, M. Vera, A. Dreszman, and L. Rey Vega, “Diffusion
assisted image reconstruction in optoacoustic tomography,” Optics
and Lasers in Engineering, vol. 178, p. 108242, 2024.
[32] K. B. Petersen and M. S. Pedersen, “The matrix cookbook,” 2008.
AP
´
ENDICE A
DEMOSTRACI
´
ON DEL LEMA 1
Mientras que (14) es una identidad est
´
andar, las identida-
des (15) y (16) requieren un an
´
alisis de factores. Es sencillo
notar que la distribuci
´
on de X|K = k es normal debido a
las hip
´
otesis del modelo. Su media y la covarianza pueden
calcularse como:
E[X|K = k] = E [E[X|Y ]|K = k] = A · µ
k
(23)
cov(X|K = k) = E [cov(X|Y )|K = k] + cov (E[X|Y ]|K = k)
= σ
2
v
· I + A · Λ
k
· A
T
(24)
De este modo, (16) puede obtenerse mediante la regla
de Bayes. Con esta t
´
ecnica es f
´
acil hallar la distribuci
´
on
conjunta:
X
Y
K=k
∼ N
Aµ
k
µ
k
,
σ
2
v
I + AΛ
k
A
T
AΛ
k
Λ
k
A
T
Λ
k
(25)
La ecuaci
´
on (15) se prueba utilizando propiedades de
variables normales multivariadas [32].
AP
´
ENDICE B
DEMOSTRACI
´
ON DEL LEMA 2
La funci
´
on score se define como ψ(x) = ∇
x
log p(x).
Por un lado, puede escribirse intercambiando el orden de
diferenciaci
´
on e integraci
´
on como:
ψ(x) =
∇
x
p(x)
p(x)
=
1
p(x)
Z
R
d
y
p(y)∇
x
p(x|y)dy (26)
En nuestro modelo X|Y = y ∼ N (Ay, σ
2
v
I), el gradiente
puede calcularse como:
∇
x
p(x|y) =
Ay − x
σ
2
v
p(x|y) (27)
y la funci
´
on score es
ψ(x) =
1
p(x)
Z
R
d
y
p(y)
Ay − x
σ
2
v
p(x|y)dy
=
A · E[Y |X = x] − x
σ
2
v
(28)
De esta manera, (17) se demuestra resolviendo (28). Por
otro lado, la densidad de probabilidad p(x) en este modelo
puede escribirse como:
p(x) =
1
(2πσ
2
v
)
d
x
2
Z
R
d
y
p(y)e
−
1
2σ
2
v
∥x−Ay∥
2
dy
=
1
(2πσ
2
v
)
d
x
2
E
e
−
1
2σ
2
v
∥x−AY ∥
2
(29)
Por lo tanto, la funci
´
on score tambi
´
en puede expresarse
como
ψ(x) = ∇
x
log E
e
−
1
2σ
2
v
∥x−AY ∥
2
(30)
Finalmente, el lema queda demostrado.
Revista elektron, Vol. 9, No. 2, pp. 76-83 (2025)
ISSN 2525-0159
83
http://elektron.fi.uba.ar
