Gaussian Splatting-协方差矩阵与分解

[toc]

2D, 3D协方差矩阵与公式6, 协方差矩阵为什么可以分解为R和S

参考:

https://stackoverflow.com/questions/65527024/create-covariance-matrix-using-ratio-and-rotation-degree

https://github.com/graphdeco-inria/gaussian-splatting/issues/392

首先, 协方差矩阵的性质:半正定,方阵, 实对称(实数+对称), 其中半正定意味着对于任何的非零向量$v$(假设$\Sigma$为$n$阶方阵, $v$为n维向量), 满足:

\[v \Sigma v^T \geq 0\]

这个公式同时也意味着协方差矩阵的所有特征值都是非负的。

举例:

\[[b_1, b_2, b_3] \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{bmatrix} [b_1, b_2, b_3]^T \geq 0\]

而协方差矩阵是对称的, 其转置等于自身, 因此$a_{ij}=a_{ji}$, 可以变为下面的形式

\[b_1^2*a_{11} + b_2^2*a_{22} + b_3^2*a_{33} + 2*b_1*b_2*a_{12} + 2*b_1*b_3*a_{13} + 2*b_2*b_3*a_{23} \geq 0\]

因为$v$是非零向量, 因此所有的$a_{ij}$都是非负的。

---

设概率空间为$(\Omega, \Sigma, P)$, 有两个随机变量序列, $X={ x_i }^{m}{i=1}$和$Y={ y_i }^{m}{i=1}$, 其中的$x_i$和$y_i$都是随机变量。 例如:

$X={x_1, x_2}$, $Y={y_1, y_2}$

其中$x_i$和$y_i$都是表示为数组的随机变量, 包含一系列的观测值:

$x_1 = [2, 4, 6], x_2 = [1, 3, 5], y_1 = [5, 10 ,15], y_2 = [3, 6, 9]$

那么对应的期望为 $E(x_i)=\int_{\Omega} x_i dP$,

协方差描述两个随机变量之间的差异分布, 协方差矩阵是两列/组随机变量之间的协方差, 是协方差的推广。

这两组随机变量之间的协方差矩阵定义为:

\[cov(X, Y) := [cov(x_i, y_j)]_{m \times n} = [ E[(x_i - u_i)(y_j - v_j)] ]_{m \times n}\] \[cov(X, Y) := [ E[(X - E(X))(Y - E(Y))^T] ]\]

协方差矩阵能够导出一个变换矩阵, 使数据完全去相关(decorrelation), 换句话说, 可以找出一组最佳的基, 以紧凑的方式表达数据。也就是主成分分析, KL-变换。参考: https://en.wikipedia.org/wiki/Covariance_matrix

---

参考: https://medium.com/@AriaLeeNotAriel/numbynum-3d-gaussian-splatting-for-real-time-radiance-field-rendering-kerbl-et-al-60c0b25e5544

渲染需要将3D Gaussian投影到2D, $J$是投影变换的仿射近似雅克比矩阵。$W$是投影变换。

\[\Sigma = JW\Sigma W^TJ^T\]

如同论文中说的那样, 最直接的方法就是创建一个对称矩阵, 创建6个变量, 然后梯度下降去优化这6个值。$a_{*}$代表其值通过对称获得。

\[\Sigma = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{*} & a_{22} & a_{23} \\ a_{*} & a_{*} & a_{33} \\ \end{bmatrix}\]

但是协方差矩阵必须是半正定的才有物理意义, 优化过程中很容易就会破坏协方差矩阵的半正定。

因此，作者提出了一种具有同样表达能力的优化表示来代替协方差矩阵, 即论文中的公式6, 介绍如何使用scale和rotation来计算协方差矩阵. 这一种各向异性的协方差表示适用于优化。公式6:

\[\Sigma=R S S^T R^T\]

为什么这么做？这种操作被称为reparametrization。利用下面的公式得到的矩阵总是满足半正定的

\[\Sigma = A^T A\]

因为对于任何的非零向量$v$, $Av$是一个维度为n的向量。

\[v^T\Sigma v = v^{T} (A^{T} A) v = (Av)^{T}Av = \|Av\|^{2} \geq 0\]

Q: 为什么$A$被进一步分解为Scaling和Rotation?

因为论文的目标就是优化每个高斯椭球的参数来得到更好的可视化效果(PSNR等指标), 而每个高斯椭球都是由均值(中心点坐标), 协方差矩阵($\Sigma$)和不透明度参数$\alpha$来定义的, 协方差矩阵控制着椭球的形状和方向. 具体来说, 是Scaling和Rotation在分别控制着椭球的形状和方向。

将其使用$RS$组合成矩阵$A$只是为了维持协方差矩阵的半正定这一物理性质, 换句话说, 使用$R$和$S$来组合成$A$其实是人为设计的策略, 我完全可以创建一个矩阵A, 把四元数$q$和缩放向量$s$的一共7个值塞进去来创建一个A; 更甚至, 我可以将法向或者颜色或者其他的什么乱七八糟的塞到$A$中, 让其参与到梯度下降中去。

下面的CUDA代码使用scale和rotation来计算协方差矩阵

mod是scale_modifier, 来自

// Forward method for converting scale and rotation properties of each
// Gaussian to a 3D covariance matrix in world space. Also takes care
// of quaternion normalization.
__device__ void computeCov3D(const glm::vec3 scale, float mod, const glm::vec4 rot, float* cov3D)
{
	// Create scaling matrix
	glm::mat3 S = glm::mat3(1.0f);  // 创建3×3零矩阵, scale赋值对角线元素
	S[0][0] = mod * scale.x;
	S[1][1] = mod * scale.y;
	S[2][2] = mod * scale.z;

	// Normalize quaternion to get valid rotation
	glm::vec4 q = rot;// / glm::length(rot);
	float r = q.x;
	float x = q.y;
	float y = q.z;
	float z = q.w;

	// Compute rotation matrix from quaternion  // 四元数转旋转矩阵
	glm::mat3 R = glm::mat3(
		1.f - 2.f * (y * y + z * z), 2.f * (x * y - r * z), 2.f * (x * z + r * y),
		2.f * (x * y + r * z), 1.f - 2.f * (x * x + z * z), 2.f * (y * z - r * x),
		2.f * (x * z - r * y), 2.f * (y * z + r * x), 1.f - 2.f * (x * x + y * y)
	);

	glm::mat3 M = S * R;  // 得到 scaling 和 rotation 矩阵的乘积

	// Compute 3D world covariance matrix Sigma
	glm::mat3 Sigma = glm::transpose(M) * M;  // M^T M 计算协方差矩阵\Sigma

	// Covariance is symmetric, only store upper right
    // 协方差矩阵对称, 只保留上三角部分
	cov3D[0] = Sigma[0][0];
	cov3D[1] = Sigma[0][1];
	cov3D[2] = Sigma[0][2];
	cov3D[3] = Sigma[1][1];
	cov3D[4] = Sigma[1][2];
	cov3D[5] = Sigma[2][2];
}

仔细观察上面的代码, 可能会有疑问: 论文中的公式6写的是:

\[\Sigma=R S S^T R^T\]

($S$的对角线元素是各个轴的缩放, 其他非对角线元素的值为0, 是对称矩阵, 所以$S=S^T$)

且附录中在介绍求导时提到:

\[M=RS\] \[\Sigma = MM^T = RSS^TR^T = RSSR^T\]

和公式6一致。

而这个代码的实现中:

\[M = SR\] \[\Sigma = M^T M = R^TS^TSR = R^TSSR\]

在 https://github.com/graphdeco-inria/gaussian-splatting/issues/762 中猜测CUDA代码中的旋转矩阵$R$与论文中的$R$并不是同一个, 是互逆的关系。 但是真的如此吗? 我们将CUDA代码中的R计算出来:

\[R = \begin{bmatrix} 1 - 2y^2 - 2z^2 & 2xy - 2rz & 2xz + 2ry \\ 2xy + 2rz & 1 - 2x^2 - 2z^2 & 2yz - 2rx \\ 2xz - 2ry & 2yz + 2rx & 1 - 2x^2 - 2y^2 \end{bmatrix} = 2 \begin{bmatrix} \frac{1}{2} - y^2 - z^2 & xy - rz & xz + ry \\ xy + rz & \frac{1}{2} - x^2 - z^2 & yz - rx \\ xz - ry & yz + rx & \frac{1}{2} - x^2 - y^2 \end{bmatrix}\]

四元数为$(q_{r}, q_{i}, q_{j}, q_{k})$, 对应代码中的$r, x, y, z$ \(R = 2 \begin{bmatrix} \frac{1}{2} - q_{j}^2 - q_{k}^2 & q_{i}q_{j} - q_{r}q_{k} & q_{i}q_{k} + q_{r}q_{j} \\ q_{i}q_{j} + q_{r}q_{k} & \frac{1}{2} - q_{i}^2 - q_{k}^2 & q_{j}q_{k} - q_{r}q_{i} \\ q_{i}q_{k} - q_{r}q_{j} & q_{j}q_{k} + q_{r}q_{i} & \frac{1}{2} - q_{i}^2 - q_{j}^2 \end{bmatrix}\)

发现与论文附录中的公式10是完全相同的

也就是说,论文中的旋转矩阵和代码中的是相同的。

真的是相同的吗？

https://github.com/graphdeco-inria/gaussian-splatting/issues/100#issuecomment-1690378854 中提到了在Camera中对旋转矩阵进行了转置, 是为了满足行主序(column-major)

https://github.com/graphdeco-inria/gaussian-splatting/issues/208 中提到了pytorch的矩阵行主序的问题

到底有什么关系？我们平时使用的opencv中的都是行主序, $M[i][j]$是第i行第j列。

让我们来用代码验证一下

在计算协方差矩阵的函数中进行输出。

__device__ void computeCov3D(const glm::vec3 scale, float mod, const glm::vec4 rot, float* cov3D)
{
	// Create scaling matrix
	glm::mat3 S = glm::mat3(1.0f);
	S[0][0] = mod * scale.x;
	S[1][1] = mod * scale.y;
	S[2][2] = mod * scale.z;

	// Normalize quaternion to get valid rotation
	glm::vec4 q = rot;// / glm::length(rot);
	float r = q.x;
	float x = q.y;
	float y = q.z;
	float z = q.w;

	// Compute rotation matrix from quaternion
	glm::mat3 R = glm::mat3(
		1.f - 2.f * (y * y + z * z), 2.f * (x * y - r * z), 2.f * (x * z + r * y),
		2.f * (x * y + r * z), 1.f - 2.f * (x * x + z * z), 2.f * (y * z - r * x),
		2.f * (x * z - r * y), 2.f * (y * z + r * x), 1.f - 2.f * (x * x + y * y)
	);
	printMat3("Rotation matrix R:", R);
	printMat3("Scaling matrix S:", S);	

	glm::mat3 M = S * R;
	printMat3("M=S*R:", M);

	glm::mat3 M2 = R * S;
	printMat3("M2=R*S:", M2);

	// Compute 3D world covariance matrix Sigma
	glm::mat3 Sigma = glm::transpose(M) * M;
	printMat3("Sigma=M^T*M:", Sigma);

	glm::mat3 Sigma2 = M * glm::transpose(M);
	printMat3("Sigma2=M*M^T:", Sigma2);

	// Covariance is symmetric, only store upper right
	cov3D[0] = Sigma[0][0];
	cov3D[1] = Sigma[0][1];
	cov3D[2] = Sigma[0][2];
	cov3D[3] = Sigma[1][1];
	cov3D[4] = Sigma[1][2];
	cov3D[5] = Sigma[2][2];
}

输出如下:

Rotation matrix R:(四元数转换得到的旋转矩阵)
0.707107 -0.707107 0.000000 
0.707107 0.707107 0.000000 
0.000000 0.000000 1.000000 

Scaling matrix S:
10.000000 0.000000 0.000000 
0.000000 20.000000 0.000000 
0.000000 0.000000 30.000000 

M=S*R:
7.071068 -14.142136 0.000000 
7.071068 14.142137 0.000000 
0.000000 0.000000 30.000000 

M2=R*S:
7.071068 -7.071068 0.000000 
14.142136 14.142137 0.000000 
0.000000 0.000000 30.000000 

Sigma=M^T*M:
250.000000 -150.000015 0.000000 
-150.000015 250.000031 0.000000 
0.000000 0.000000 900.000000 

Sigma2=M*M^T:
100.000008 0.000001 0.000000 
0.000001 400.000031 0.000000 
0.000000 0.000000 900.000000

可以看到:

1. $SR$与$RS$的结果并不同, 旋转与缩放的顺序是有影响的
2. glm的是列主序, 先列后行。即我们读取的$mat[i][j]$读取的是第i列第j行, 因此代码中实现的$M = S * R$, 对应到论文中, 其实是$M=S^TR^T$, 也就是$(RS)^T$, 最终的协方差矩阵就是$M^TM=RS(RS)^T=RSS^TR^T=RSSR^T$, 与论文中的公式相同。

---

旋转和缩放的顺序问题

先说结论:

这个问题其实取决于你的情况:

1. 如果你的缩放只是一个标量(或者说你在空间中所有轴上的缩放程度相同), 那么先缩放再旋转和先旋转再缩放是一样的。
2. 但是如果三个轴的缩放不完全相同, 那么两个操作的顺序就该被慎重考虑。

旋转和缩放的顺序是否会影响结果?$RS是否等于SR$? 在回答这个问题以前, 很容易忽略的一个问题是, 你所说的”缩放”(scaling), 在空间的各个轴上的缩放程度是否相同? 换句话说: scaling矩阵$S$的对角线元素是否是完全相同的? 这影响着$R$和$S$是否可以交换顺序.

\[S = \begin{bmatrix} s_1 & 0 & 0 \\ 0 & s_2 & 0 \\ 0 & 0 & s_3 \end{bmatrix}\]

先旋转 (R), 后缩放 (S)

一个向量 (v) 先旋转后缩放，然后 (S) 对旋转后的向量进行缩放。 [ v’ = SRv ] 这表示先通过 (R) 转换得到的矩阵(v)，然后在这个新的状态下对这个矩阵进行 (S) 缩放。

先缩放 (S), 后旋转 (R)

(v) 先缩放，然后 (R) 对缩放后的向量进行旋转。 [ v’’ = RSv ] 这表示先在原始矩阵状态下进行缩放矩阵 (S) 缩放，然后 (R) 对缩放后的向量进行旋转。

---

Q:$RSSR^T$与$R^TSSR$的值相同吗?

依然是上面的测试代码, 这次将scale的3个值都设置为10.0f, 计算的结果如下

Sigma=M^T*M:
100.000008 -0.000000 0.000000 
-0.000000 100.000008 0.000000 
0.000000 0.000000 100.000000 
Sigma2=M*M^T:
100.000008 0.000000 0.000000 
0.000000 100.000008 0.000000 
0.000000 0.000000 100.000000 

可以发现, $RSSR^T$与$R^TSSR$的值相同, 从物理意义上也能解释的通, scale控制的其实是椭球在各个方向的半径, 当scale相同时, 此时是一个标准的圆球, 旋转也就失去了意义。

结论: vec3参数scale的$s_1, s_2, s_3$不是相同的, 那么$RSSR^T$与$R^TSSR$结果不同; 但是, 当对角线元素相同, 甚至没有进行缩放, 那么此时$RSSR^T$与$R^TSSR$的值就是相同的。