如何在黎曼意义下定义相关矩阵的内均值？
如何在黎曼意义下定义相关矩阵的内均值？
正文
在今天的这篇推文中，只有问题，没有答案。
PSD锥（协方差矩阵的集合）的黎曼几何形状非常好理解，大家可以参考下面的两个课件：
https://www.covariance2017.eu/storage/Part_I_Talk_ICCV_2017_v3.pdf
http://indico.ictp.it/event/a08167/session/124/contribution/85/material/0/0.pdf
通常考虑的内积（平滑取决于点）为：
除了具有许多不变性外，该黎曼矩阵还具有一些特殊的特征，即它的关联距离等于Fisher-Rao距离（在高斯概率密度空间上的黎曼度量）的两倍。
这两个黎曼矩阵在非常相似的空间（高斯分布的密度是用相同的均值进行参数化，协方差可以识别到协方差矩阵中）上的这种联系通过Fréchet-Darmois-Cramér-Rao不等式给出了很好的统计解释：
黎曼矩阵引起的协方差矩阵空间的曲率是统计估计不确定性的简单函数。
确实，Fréchet-Darmois-Cramér-Rao不等式从本质上说，对于一个无偏估计量，其方差受Fisher信息倒数的下限（Fisher-Rao Riemannian矩阵是Fisher信息矩阵的二次形式）：Fisher信息越高，方差越小（反之）。这意味着，对于难以估计的值（高方差），空间相当平坦；而对于易于估计的值（低方差），空间更弯曲。
为了更精确地说明这个想法，让我们们考虑以下小示例：密度以2 x 2相关矩阵参数化的二元中心高斯分布。
首先，得出Fisher信息：
从Fisher信息中，我们可以得到相关估计量方差的Fréchet–Darmois–Cramér–Rao下界：
我们们显示以下值。绝对相关性越高，估算方差的下限越低。
import operator
import numpy as np
from scipy.linalg import sqrtm
from scipy.linalg import fractional_matrix_power
from numpy.linalg import inv, eig
from matplotlib import pyplot as plt
from matplotlib import animation
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from pprint import pprint
def fisher_information(rho):
return ((rho - 1)**2 * (rho + 1)**2) / (3 * (rho**2 + 1))
rhos = np.linspace(-1, 1, 1000)
plt.plot(rhos, [fisher_information(rho) for rho in rhos])
plt.title(r'$var(\widehat{\rho}) \geq ' +
r'\dfrac{(\rho - 1)^2(\rho + 1)^2}{3(\rho^2 + 1)}$',
fontsize=16)
plt.show()
备注：估计低（绝对）值的相关性时，不确定性非常高：。取中值的系数的标准偏差很大！
现在，我们将可视化表示两个相关矩阵之间的距离。
我们在任意两个相关矩阵之间的所有成对距离的表面在下方显示。
# Riemannian distance(A, B) = (sum_i ln(lambda_i)^2)^(1/2),
# where lambda_i are the A^(-1) * B
def riemann_dist(A, B):
eigenvals, eigenvecs = eig(inv(A).dot(B))
return np.sqrt(np.sum(np.log(eigenvals)**2))
N = 500
Z = np.zeros((N, N))
for i, rho_1 in enumerate(np.linspace(-0.99, 0.99, N)):
corr_1 = np.array([[1, rho_1],
[rho_1, 1]])
for j, rho_2 in enumerate(np.linspace(-0.99, 0.99, N)):
corr_2 = np.array([[1, rho_2],
[rho_2, 1]])
dist = 0.5 * riemann_dist(corr_1, corr_2)
Z[i, j] = dist
X = np.outer(np.linspace(-0.99, 0.99, N), np.ones(N))
Y = X.copy().T
fig = plt.figure(figsize=(15, 10))
ax = plt.axes(projection='3d')
ax.plot_surface(X, Y, Z, rcount=100, ccount=100,
cmap='viridis', edgecolor='none',
linewidth=0, antialiased=True, alpha=0.9)
ax.set_title(r'Pairwise Riemannian distances between ' +
r'$2 \times 2$ correlation matrices' + "\n",
fontsize=16)
plt.xlabel(r'$\rho_1$', fontsize=16)
plt.ylabel(r'$\rho_2$', fontsize=16)
plt.show()
对于高（绝对）相关值，与应用于低（绝对）相关值的相同小变化相比，相关值中的小变化会导致距离的较大变化。在较高（绝对）相关值时，空间更加弯曲。
为什么统计学家会更喜欢这种黎曼矩阵：与平面空间（例如欧几里得或Wasserstein）相比，空间曲率可以更好地区分信号和噪声（由于统计估计）。
def init():
ax.plot_surface(X, Y, Z, rcount=50, ccount=50,
cmap='viridis', edgecolor='none')
return fig,
def animate(i):
ax.view_init(elev=10, azim=i*4)
return fig,
anim = animation.FuncAnimation(
fig, animate, init_func=init, frames=90, interval=50,
blit=True)
fn = 'pairwise_corr_dists'
anim.save(fn + '.mp4', writer='ffmpeg', fps=1000 / 50)
anim.save(fn + '.gif', writer='imagemagick', fps=1000 / 50)
plt.rcParams['animation.html'] = 'html5'
anim
在讲述了黎曼矩阵的使用并讨论了其统计解释之后，回到最初的问题：如何定义相关矩阵的内在黎曼均值？
我们阅读了几篇电气工程论文，这些论文在处理相关矩阵（此集合有时称为椭圆）时，使用了以下事实：它们是协方差矩阵的子集（此集合也称为正半定（PSD）锥），然后利用PSD锥的几何形状及其黎曼矩阵得出均值、中位数，以及从相关矩阵中得出其他几何量。
但是，我们不清楚这种方法是否有效。在下面的二维案例中，我们将说明为什么我们认为这不一定是最好的方法。对于搞数学的人来说，这可以简明表达：
相关矩阵的子流形（其中由引起的黎曼矩阵）不是完全测地子流形，即中的测地线不一定是中的测地线。
在下面的动画中可以很容易地看到这一点：
对于形状为的PSD矩阵，相关矩阵（椭圆）被限制为一个简单的分段（x = 1，y = 1，z = -1..1）（以橙色显示）。
让我们们考虑和两个相关矩阵。当约束在椭圆（橙色段）上时，和之间的测地线是和之间的子段。
但是，当将和作为中的点（即协方差矩阵）时，和之间的测地线是绿色曲线。
因此，并不完全是测地线。
关于均值。两个相关矩阵的黎曼均值是测地线（）的中点（或，其中是黎曼距离，即一般Fréchet均值定义计算超过两个点的均值），并在下面显示为绿色点。两个相关矩阵的均值通常不是相关矩阵，而是协方差矩阵。
如果我们们仅希望或需要使用相关矩阵，该怎么办？
论文通常通过其方差对平均协方差进行归一化，以获得均值相关性，即，由下面的绿色三角形显示。
将平均协方差投影到相关空间的一种更几何的方法是找到相对于该平均协方差的黎曼距离d最接近的相关矩阵，即，这里，。该最接近的相关矩阵在下面显示为红色三角形。
寻找的相关矩阵解。它在下面显示为洋红色点，并且测地线从该点到（洋红色）和（红色）。
我们相信2.和3.是等效的。
请注意，通常，方法1.和2.(或3.)不会产生相同的“均值”相关矩阵。
问题：
- “黎曼”平均相关矩阵到底应该是什么？我们倾向于2.或3.。
- 一个定义是否提供更好的属性？
- 这些属性是什么？
- 我们们是否可以定义测地线停留在椭圆上的内在黎曼均值？（不是3.）
下面的动画很好地说明了问题：
def init():
ax.plot_wireframe(X, Y, Z, ccount=20, rcount=20,
label='boundaries of the PSD cone ' +
'(set of covariance matrices)')
return fig,
fig = plt.figure(figsize=(15, 10))
ax = Axes3D(fig)
# SPD cone
x = np.linspace(0, 2, num=500)
y = np.linspace(0, 2, num=500)
X, Y = np.meshgrid(x, y)
Z = np.sqrt(X * Y)
ax.plot_wireframe(X, Y, Z, ccount=20, rcount=20,
label='boundaries of the PSD cone ' +
'(set of covariance matrices)')
ax.plot_wireframe(X, Y, -Z, ccount=20, rcount=20)
# Elliptope
zline = np.linspace(-1, 1, 1000)
xline = [1] * len(zline)
yline = [1] * len(zline)
ax.plot3D(xline, yline, zline, 'orange',
label='Elliptope (subset of correlation matrices)')
A = np.array([[1, -0.4],
[-0.4, 1]])
B = np.array([[1, 0.95],
[0.95, 1]])
ax.scatter(A[0,0],
A[1,1],
A[0,1],
c='b', marker='o', s=30,
label='A (x=1, y=1, z=-0.4)')
ax.scatter(B[0,0],
B[1,1],
B[0,1],
c='k', marker='o', s=30,
label='B (x=1, y=1, z=0.95)')
# Geodesic between A and B
# p(t) = A^(1/2) * (A^(-1/2) * B * A^(-1/2))^t * A^(1/2)
sqrt_A = fractional_matrix_power(A, 0.5)
inv_sqrt_A = fractional_matrix_power(A, -0.5)
x_geod = []
y_geod = []
z_geod = []
for t in np.linspace(0, 1, 100):
geod_point = (
sqrt_A.dot(fractional_matrix_power(
inv_sqrt_A.dot(B.dot(inv_sqrt_A)), t)
.dot(sqrt_A)))
x_geod.append(geod_point[0, 0])
y_geod.append(geod_point[1, 1])
z_geod.append(geod_point[0, 1])
ax.plot3D(x_geod, y_geod, z_geod, 'green',
label=r'Geodesic $\gamma(t), t \in [0, 1],' +
'\gamma(0) = A, \gamma(1) = B$')
# Riemannian mean
# the mean of two correlation matrix is not a correlation matrix in general
t = 0.5
riemann_mean = (sqrt_A.dot(fractional_matrix_power(
inv_sqrt_A.dot(B.dot(inv_sqrt_A)), t).dot(sqrt_A)))
ax.scatter(riemann_mean[0,0],
riemann_mean[1,1],
riemann_mean[0,1],
c='g', marker='o', s=30,
label=r'Riemannian mean(A, B) = $\gamma(0.5)$')
# Riemannian mean, which is a covariance,
# "projected" back to the correlation space
# diag(Sigma)^(-1/2) * Sigma * diag(Sigma)^(-1/2)
corr_cvt = np.diag(np.sqrt(np.diag(riemann_mean)**(-1))).dot(
riemann_mean.dot(
np.diag(np.sqrt(np.diag(riemann_mean)**(-1)))))
ax.scatter(corr_cvt[0,0],
corr_cvt[1,1],
corr_cvt[0,1],
c='g', marker='^', s=30,
label='Riemannian mean(A, B) normalized by variance')
# point in the correlation space which minimizes the distance
# (in the Riemannian metric sense for the SPD space) to the covariance mean
closest_corr = min([(riemann_dist(np.array([[1, t], [t, 1]]),
riemann_mean), np.array([[1, t], [t, 1]]))
for t in np.linspace(-0.99, 0.99, 5000)],
key=operator.itemgetter(0))[1]
ax.scatter(closest_corr[0,0],
closest_corr[1,1],
closest_corr[0,1],
c='r', marker='^', s=30,
label='Closest correlation matrix to the Riemannian mean(A, B)'
+ '\n' + 'wrt the Riemannian distance')
# Fréchet Mean:
# m = argmin_{C \in E} \sum_{i=1}^N d^2(C, C_i).
frechet_mean = min(
[(riemann_dist(np.array([[1, t], [t, 1]]), A)**2 +
riemann_dist(np.array([[1, t], [t, 1]]), B)**2,
np.array([[1, t], [t, 1]]))
for t in np.linspace(-0.99, 0.99, 5000)],
key=operator.itemgetter(0))[1]
ax.scatter(frechet_mean[0,0],
frechet_mean[1,1],
frechet_mean[0,1],
c='magenta', marker='o', s=30,
label='Fréchet Mean M(A, B)')
# Geodesic between A and Frechet Mean M
# p(t) = A^(1/2) * (A^(-1/2) * M * A^(-1/2))^t * A^(1/2)
sqrt_A = fractional_matrix_power(A, 0.5)
inv_sqrt_A = fractional_matrix_power(A, -0.5)
x_geod = []
y_geod = []
z_geod = []
for t in np.linspace(0, 1, 100):
geod_point = (
sqrt_A.dot(fractional_matrix_power(
inv_sqrt_A.dot(frechet_mean.dot(inv_sqrt_A)), t)
.dot(sqrt_A)))
x_geod.append(geod_point[0, 0])
y_geod.append(geod_point[1, 1])
z_geod.append(geod_point[0, 1])
ax.plot3D(x_geod, y_geod, z_geod, 'magenta',
label=r'Geodesic $\gamma(t), t \in [0, 1],' +
'\gamma(0) = A, \gamma(1) = M$')
# Geodesic between B and Frechet Mean M
# p(t) = B^(1/2) * (B^(-1/2) * M * B^(-1/2))^t * B^(1/2)
sqrt_A = fractional_matrix_power(B, 0.5)
inv_sqrt_A = fractional_matrix_power(B, -0.5)
x_geod = []
y_geod = []
z_geod = []
for t in np.linspace(0, 1, 100):
geod_point = (
sqrt_A.dot(fractional_matrix_power(
inv_sqrt_A.dot(frechet_mean.dot(inv_sqrt_A)), t)
.dot(sqrt_A)))
x_geod.append(geod_point[0, 0])
y_geod.append(geod_point[1, 1])
z_geod.append(geod_point[0, 1])
ax.plot3D(x_geod, y_geod, z_geod, 'red',
label=r'Geodesic $\gamma(t), t \in [0, 1],' +
'\gamma(0) = B, \gamma(1) = M$')
# Euclidean mean
eucl_mean = (A + B) / 2
ax.scatter(eucl_mean[0,0],
eucl_mean[1,1],
eucl_mean[0,1],
c='orange', marker='o', s=30,
label='Euclidean mean(A, B)')
ax.legend()
anim = animation.FuncAnimation(
fig, animate, init_func=init, frames=90, interval=50,
blit=True)
fn = 'SPD_cone'
anim.save(fn + '.mp4', writer='ffmpeg', fps=1000 / 50)
anim.save(fn + '.gif', writer='imagemagick', fps=1000 / 50)
plt.rcParams['animation.html'] = 'html5'
anim
