45
45
#TODO:
46
46
# - Add some documentation to `proximal`, `gradient` and `convex_conj`
47
47
# (use See Also when applicable, otherwise a full doc)
48
- # - Unify citations
49
48
50
49
51
50
class LpNorm (Functional ):
@@ -1222,7 +1221,7 @@ def convex_conj(self):
1222
1221
1223
1222
@property
1224
1223
def proximal (self ):
1225
- """Return the proximal factory of the functional.
1224
+ """A proximal factory for this functional.
1226
1225
1227
1226
This is the zero operator.
1228
1227
"""
@@ -1256,37 +1255,44 @@ def __repr__(self):
1256
1255
allow_mixed_seps = False )
1257
1256
1258
1257
1259
- #TODO: continue here
1260
-
1261
1258
class KullbackLeibler (Functional ):
1262
1259
1263
1260
r"""The Kullback-Leibler divergence functional.
1264
1261
1265
1262
Notes
1266
1263
-----
1267
- The functional :math:`F` with prior :math:`g>=0` is given by:
1264
+ The Kullback-Leibler divergence with prior :math:`g>=0` is defined as
1268
1265
1269
1266
.. math::
1270
- F (x)
1271
- =
1267
+ \text{KL} (x)
1268
+ & =
1272
1269
\begin{cases}
1273
- \sum_{i} \left( x_i - g_i + g_i \log \left( \frac{g_i}{x_i}
1274
- \right) \right) & \text{if } x_i > 0 \forall i
1270
+ \sum_{i} \left( x_i - g_i + g_i \ln \left( \frac{g_i}{x_i}
1271
+ \right) \right) & \text{if } x_i > 0 \text{ for all } i,
1275
1272
\\
1276
- +\infty & \text{else .}
1273
+ +\infty & \text{otherwise .}
1277
1274
\end{cases}
1275
+ \quad (\mathbb{R}^n\text{-like space}) \\[2ex]
1276
+ \text{KL}(x)
1277
+ &=
1278
+ \begin{cases}
1279
+ \int \left(
1280
+ x(t) - g(t) + g(t) \ln\left(\frac{g(t)}{x(t)}\right)
1281
+ \right)\, \mathrm{d}t & \text{if } x(t) > 0 \text{ for all } t,
1282
+ \\
1283
+ +\infty & \text{otherwise.}
1284
+ \end{cases}
1285
+ \quad (L^p-\text{like space})
1278
1286
1279
- Note that we use the common definition 0 log(0) := 0.
1280
- KL based objectives are common in MLEM optimization problems and are often
1287
+ Note that we use the common convention :math:`0 \ln 0 := 0` .
1288
+ KL- based objectives are common in MLEM optimization problems and are often
1281
1289
used as data-matching term when data noise governed by a multivariate
1282
1290
Poisson probability distribution is significant.
1283
1291
1284
- The functional is related to the Kullback-Leibler cross entropy functional
1285
- `KullbackLeiblerCrossEntropy`. The KL cross entropy is the one
1286
- diescribed in `this Wikipedia article
1287
- <https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence>`_, and
1288
- the functional :math:`F` is obtained by switching place of the prior and
1289
- the varialbe in the KL cross entropy functional.
1292
+ This functional is related to the `KullbackLeiblerCrossEntropy`
1293
+ described in `this Wikipedia article
1294
+ <https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence>`_,
1295
+ in that they have flipped roles of variable :math:`x` and prior :math:`g`.
1290
1296
1291
1297
For a theoretical exposition, see `[Csiszar1991]
1292
1298
<http://www.jstor.org/stable/2241918>`_.
@@ -1340,23 +1346,23 @@ def prior(self):
1340
1346
def _call (self , x ):
1341
1347
"""Return ``self(x)``.
1342
1348
1343
- If any components of ``x`` is non-positive, the value is positive
1349
+ If any component of ``x`` is non-positive, the value is positive
1344
1350
infinity.
1345
1351
"""
1346
1352
# Lazy import to improve `import odl` time
1347
1353
import scipy .special
1348
1354
1349
1355
if self .prior is None :
1350
- tmp = (( x - 1 - np .log (x )).inner (self .domain .one () ))
1356
+ integral = (x - 1 - np .log (x )).inner (self .domain .one ())
1351
1357
else :
1352
- tmp = ( (x - self .prior +
1353
- scipy .special .xlogy (self .prior , self .prior / x ))
1354
- .inner (self .domain .one () ))
1355
- if np .isnan (tmp ):
1358
+ integrand = (x - self .prior +
1359
+ scipy .special .xlogy (self .prior , self .prior / x ))
1360
+ integral = integrand .inner (self .domain .one ())
1361
+ if np .isnan (integral ):
1356
1362
# In this case, some element was less than or equal to zero
1357
1363
return np .inf
1358
1364
else :
1359
- return tmp
1365
+ return integral
1360
1366
1361
1367
@property
1362
1368
def gradient (self ):
@@ -1365,7 +1371,7 @@ def gradient(self):
1365
1371
For a prior :math:`g` is given by
1366
1372
1367
1373
.. math::
1368
- \nabla F (x) = 1 - \frac{g}{x}.
1374
+ \nabla \text{KL} (x) = 1 - \frac{g}{x}.
1369
1375
1370
1376
The gradient is not defined if any component of :math:`x` is
1371
1377
non-positive.
@@ -1392,7 +1398,7 @@ def _call(self, x):
1392
1398
1393
1399
@property
1394
1400
def proximal (self ):
1395
- """Return the `proximal factory` of the functional.
1401
+ """A `proximal factory` for this functional.
1396
1402
1397
1403
See Also
1398
1404
--------
@@ -1401,12 +1407,17 @@ def proximal(self):
1401
1407
odl.solvers.nonsmooth.proximal_operators.proximal_convex_conj :
1402
1408
Proximal of the convex conjugate of a functional.
1403
1409
"""
1404
- return proximal_convex_conj (proximal_convex_conj_kl ( space = self . domain ,
1405
- g = self .prior ))
1410
+ return proximal_convex_conj (
1411
+ proximal_convex_conj_kl ( space = self . domain , g = self .prior ))
1406
1412
1407
1413
@property
1408
1414
def convex_conj (self ):
1409
- """The convex conjugate functional of the KL-functional."""
1415
+ """The convex conjugate of the KL functional.
1416
+
1417
+ See Also
1418
+ --------
1419
+ KullbackLeiblerConvexConj
1420
+ """
1410
1421
return KullbackLeiblerConvexConj (self .domain , self .prior )
1411
1422
1412
1423
def __repr__ (self ):
@@ -1429,24 +1440,38 @@ def __repr__(self):
1429
1440
1430
1441
class KullbackLeiblerConvexConj (Functional ):
1431
1442
1432
- r"""The convex conjugate of Kullback-Leibler divergence functional.
1443
+ r"""The convex conjugate of the Kullback-Leibler divergence functional.
1433
1444
1434
1445
Notes
1435
1446
-----
1436
- The functional :math:`F^*` with prior :math:`g > 0` is given by
1447
+ The convex conjugate :math:`\text{KL}^*` of the KL divergence is given
1448
+ by
1449
+
1450
+ The Kullback-Leibler divergence with prior :math:`g>=0` is defined as
1437
1451
1438
1452
.. math::
1439
- F^*(x) =
1453
+ \text{KL}(x)
1454
+ &=
1440
1455
\begin{cases}
1441
- \sum_{i} \left( -g_i \ln(1 - x_i) \right)
1442
- & \text{if } x_i < 1 \forall i
1456
+ \sum_{i} \left( -g_i \ln(1 - x_i) \right) & \text{if }
1457
+ x_i < 1 \text{ for all } i,
1443
1458
\\
1444
- +\infty & \text{else }
1459
+ +\infty & \text{otherwise. }
1445
1460
\end{cases}
1461
+ \quad (\mathbb{R}^n\text{-like space}) \\[2ex]
1462
+ \text{KL}(x)
1463
+ &=
1464
+ \begin{cases}
1465
+ \int \big(-g(t)\ln\left(1 - x(t)\big)
1466
+ \right)\, \mathrm{d}t & \text{if } x(t) < 1 \text{ for all } t,
1467
+ \\
1468
+ +\infty & \text{otherwise.}
1469
+ \end{cases}
1470
+ \quad (L^p-\text{like space})
1446
1471
1447
1472
See Also
1448
1473
--------
1449
- KullbackLeibler : convex conjugate functional
1474
+ KullbackLeibler : convex conjugate
1450
1475
"""
1451
1476
1452
1477
def __init__ (self , space , prior = None ):
@@ -1472,37 +1497,42 @@ def __init__(self, space, prior=None):
1472
1497
1473
1498
@property
1474
1499
def prior (self ):
1475
- """The prior in convex conjugate Kullback-Leibler functional."""
1500
+ """The prior in the convex conjugate of the KL functional."""
1476
1501
return self .__prior
1477
1502
1478
1503
# TODO(#440): use integration operator when available
1479
1504
def _call (self , x ):
1480
1505
"""Return ``self(x)``.
1481
1506
1482
- If any components of ``x`` is larger than or equal to 1, the value is
1507
+ If any component of ``x`` is larger than or equal to 1, the value is
1483
1508
positive infinity.
1484
1509
"""
1485
1510
# Lazy import to improve `import odl` time
1486
1511
import scipy .special
1487
1512
1488
1513
if self .prior is None :
1489
- tmp = self .domain .element (
1490
- - 1.0 * (np .log (1 - x ))).inner (self .domain .one ())
1514
+ integral = (- 1.0 * (np .log1p (- x ))).inner (self .domain .one ())
1491
1515
else :
1492
- tmp = self . domain . element ( - scipy .special .xlogy (
1493
- self . prior , 1 - x )) .inner (self .domain .one ())
1494
- if np .isnan (tmp ):
1516
+ integrand = - scipy .special .xlog1py ( self . prior , - x )
1517
+ integral = integrand .inner (self .domain .one ())
1518
+ if np .isnan (integral ):
1495
1519
# In this case, some element was larger than or equal to one
1496
1520
return np .inf
1497
1521
else :
1498
- return tmp
1522
+ return integral
1499
1523
1500
1524
@property
1501
1525
def gradient (self ):
1502
- """Gradient operator of the functional.
1526
+ """Gradient operator of this functional.
1503
1527
1504
- The gradient is not defined in points where one or more components
1505
- are larger than or equal to one.
1528
+ The gradient of the convex conjugate of the KL divergence is given
1529
+ by
1530
+
1531
+ .. math::
1532
+ \n abla \t ext{KL}^*(x) = \f rac{g}{1 - x}.
1533
+
1534
+ The gradient is not defined in points where any component of :math:`x`
1535
+ is (larger than or) equal to one.
1506
1536
"""
1507
1537
functional = self
1508
1538
@@ -1530,7 +1560,7 @@ def _call(self, x):
1530
1560
1531
1561
@property
1532
1562
def proximal (self ):
1533
- """Return the `proximal factory` of the functional.
1563
+ """A `proximal factory` for this functional.
1534
1564
1535
1565
See Also
1536
1566
--------
@@ -1543,7 +1573,10 @@ def proximal(self):
1543
1573
1544
1574
@property
1545
1575
def convex_conj (self ):
1546
- """The convex conjugate functional of the conjugate KL-functional."""
1576
+ """The convex conjugate functional of the KL convex conjugate.
1577
+
1578
+ This is the original KL divergence.
1579
+ """
1547
1580
return KullbackLeibler (self .domain , self .prior )
1548
1581
1549
1582
def __repr__ (self ):
@@ -1556,6 +1589,7 @@ def __repr__(self):
1556
1589
allow_mixed_seps = False )
1557
1590
1558
1591
1592
+ #TODO: continue here
1559
1593
class KullbackLeiblerCrossEntropy (Functional ):
1560
1594
1561
1595
r"""The Kullback-Leibler Cross Entropy divergence functional.
@@ -1568,7 +1602,7 @@ class KullbackLeiblerCrossEntropy(Functional):
1568
1602
F(x)
1569
1603
=
1570
1604
\begin{cases}
1571
- \sum_{i} \left( g_i - x_i + x_i \log \left( \frac{x_i}{g_i}
1605
+ \sum_{i} \left( g_i - x_i + x_i \ln \left( \frac{x_i}{g_i}
1572
1606
\right) \right)
1573
1607
& \text{if } g_i > 0 \forall i
1574
1608
\\
@@ -1638,7 +1672,7 @@ def prior(self):
1638
1672
def _call (self , x ):
1639
1673
"""Return ``self(x)``.
1640
1674
1641
- If any components of ``x`` is non-positive, the value is positive
1675
+ If any component of ``x`` is non-positive, the value is positive
1642
1676
infinity.
1643
1677
"""
1644
1678
# Lazy import to improve `import odl` time
@@ -2204,7 +2238,7 @@ def __repr__(self):
2204
2238
2205
2239
class NuclearNorm (Functional ):
2206
2240
2207
- r"""Nuclear norm for matrix valued functions.
2241
+ r"""Nuclear norm for matrix- valued functions.
2208
2242
2209
2243
Notes
2210
2244
-----
@@ -2294,8 +2328,7 @@ def _asvector(self, arr):
2294
2328
2295
2329
def _call (self , x ):
2296
2330
"""Return ``self(x)``."""
2297
-
2298
- # Convert to array with most
2331
+ # Convert to array with "outer" indices last
2299
2332
arr = self ._asarray (x )
2300
2333
svd_diag = np .linalg .svd (arr , compute_uv = False )
2301
2334
0 commit comments