雅可比矩阵及应用

2018年6月21日 | Oscar | OpenGL , 实时渲染 , 数学

雅可比矩阵及应用

雅可比矩阵

　　假设某函数 $\mathbb{R}^n$ 需要映射到另一个空间 $\mathbb{R}^m$ 中，雅克比矩阵就是从 $\mathbb{R}^n$ 到 $\mathbb{R}^m$ 的线性映射，其重要意义在于它表现出了多维对多维空间的一个最佳线性估计。因此，雅可比矩阵类似于单变量函数中的导数。事实上，在单变量函数中，导数就是 $1\times 1$ 阶的雅可比矩阵。
　　注意，以下的推导的矩阵都是行为主（row major），与 OpenGL 不同，而与 DX 相同。由于代码示例使用的是 OpenGL(GLSL)，所以代码中使用的是列为主（colomn major）。
　

雅可比矩阵的直观感觉

　　对于一个单变量的函数，如 $f:\mathbb R \to \mathbb R$ ，有

$f^\prime(x) = \lim_{\Delta x \to 0}\frac{f(x + \Delta x) – f(x)}{\Delta x}.$

　　根据极限思想，就有
$f(x + \Delta x) \approx f(x) + f'(x) \Delta x$
　　上述就是基本的 1×1 维导数的思想，即用 $x$ 的一维导数来线性估计 $f(x)$ 空间的变化情况。
　　如果将 $f$ 扩展到 $m\times n$ 维的情况， $f:\mathbb R^n \to \mathbb R^m$ ：
$f(\underbrace{x}_ {n \times 1}+\underbrace{\Delta x}_ {n \times 1}) \approx\underbrace{f(x)}_ {m \times 1}+\underbrace{f^\prime(x)}_ {?}\underbrace{\Delta x}_ {n\times 1}$
　　如果尝试用 $\mathbb R^n$ 来线性估计 $\mathbb R^m$ 的值时，就需要 $f^\prime(x)$ 表现为 $m\times n$ 的矩阵，注意 $f(x)$ 是 $m \times 1$ 维，而 $x$ 是 $n \times 1$ 维的。在这种思想下，自然而然会将多维的偏微分拓展为 $m\times n$ 的矩阵，将所有的偏微分以 row major 的形式排列成矩阵，也就成了雅可比矩阵。
　　实际上，雅可比矩阵只是偏导数的集合，本质上与偏导数相差无几。单一的偏导数是对单一维度的考量，可以得到解析式的情况下，计算 $x$ 偏导数，往往会得到包含 $(y,z\dots)$ 的方程式；而雅可比行列式的运算时各个方向偏导数的叠加，根据输入方向的不同，可以得到线性估计值。
　　雅克比矩阵相乘的拆项形式，本质上也只是偏导数后相加，概念并不难理解。
　

雅可比矩阵的一个应用

　　在渲染管线进行阴影渲染的时候，会同时向片元着色器推送两个空间的顶点信息。

相机本次 draw call 的渲染顶点数据
手动在顶点着色器转换到阴影（光照空间）的顶点数据以及提前渲染好的阴影贴图。

　
　　在一些比较高级的阴影算法中，比如 PCSS或者 CSM 等，需要判断着色点邻近位置的深度情况。这就是一个 $2\times 2$ 维空间映射问题。
　　设本帧渲染的图片的二维坐标轴为 $(x,y)$ ，光照空间渲染图片的二维坐标轴为 $(u,v)$ 。需要判断 $(x,y)$ 空间中相对光源垂直方向上一些样本点的阴影状况，而垂直方向的深度信息是存储在 $(u,v)$ 空间中的。
　　但是重点不是单纯点对点的空间映射，需要根据 $z_{uv}$ （后面记为 $z$ ，注意，本次 draw call 的深度值始终没有出现过。）相对 $(x,y)$ 空间的偏微分，来线性估计 $z_{uv}$ 的变化情况，从而调整 $(x,y)$ 坐标，获得 $(x+\Delta x,y+\Delta y)$ 在光照空间的偏移深度值。
　　实际上，如果不考虑性能损耗，多渲染管线，偏移一次进行一次进行偏微分操作，得到的结果会比线性估计值更好。不过，采样 10 个点，需要进行 10 次渲染管线的顶点操作，这是不可接受的。这么做也是没有意义的浪费。
　　在 glsl 中，已有 dFdx 和 dFdy 来进行偏微分（梯度）操作，可以得到 $\partial u\over \partial x$ 、 $\partial u\over \partial y$ 、 $\partial u\over \partial z$ 、 $\partial v\over \partial x$ 、 $\partial v\over \partial y$ 和 $\partial v\over \partial z$ ，需要求的则是 $\partial z\over \partial u$ 和 $\partial z\over \partial v$ ，也就是 $(x,y)$ 空间中，深度相对于 $(u,v)$ 空间的偏移量。
　
　　首先，是阴影（光照）空间的手动透视除法工作，获得 $(u,v)$ 坐标。

  vec2 uv = LightSpacePos.xy / LightSpacePos.w; 
  float z = LightSpacePos.z / LightSpacePos.w;// 深度

　　然后，进行偏微分操作，得到上述所有已知变量。

vec3 duvz_dx = dFdx(vec3(uv,z));
vec3 duvz_dy = dFdy(vec3(uv,z));

　　这个时候可以构建雅可比矩阵：
$\begin{bmatrix} {\partial u\over \partial x}&{\partial u\over \partial y} \ {\partial v\over \partial x}&{\partial v\over \partial y} \end{bmatrix}$
　　事实上，这里也是最基础的线性变换。以 $(x,y)$ 为基空间的矩阵，是将点从其他空间转换到基空间中，所以这里需要进行逆矩阵操作。即：
${ \begin{bmatrix} {\partial z\over \partial u}&{\partial z\over \partial v} \end{bmatrix} }={ \begin{bmatrix} {\partial z\over \partial u}&{\partial z\over \partial v} \end{bmatrix} } \cdot { \begin{bmatrix} {\partial u\over \partial x}&{\partial u\over \partial y} \ {\partial v\over \partial x}&{\partial v\over \partial y} \end{bmatrix} }^{-1}$
　　二阶的逆矩阵并不难求，可以手动写在代码里。可以手动对偏微分进行运算，它是符合运算规则的。

vec2 dz_duv = vec2(0.0, 0.0);

dz_duv.x = duvz_dy.y * duvz_dx.z;
dz_duv.x -= duvz_dx.y * duvz_dy.z;

dz_duv.y = duvz_dx.x * duvz_dy.z;
dz_duv.y -= duvz_dy.x * duvz_dx.z;

float det = (duvz_dx.x * duvz_dy.y) - (duvz_dx.y * duvz_dy.x);

dz_duv /= det;

　　在使用的时候，给定基准值 $z_0 = z_{uv}$ ，就可以根据样本偏移长度来估计该点在 $(u,v)$ 空间的深度值。这个样本偏移值的采用通常会使用泊松分布。

float BiasedZ(
  float z0,
  vec2 dz_duv,
  vec2 offset)
{
  return z0 + dot(dz_duv, offset);
}

雅可比行列式

　　由上面可以看出，雅可比矩阵也只是一堆偏导数的集合，性质表现为一个更全面的导数，之所以雅可比矩阵会被总结起来，也与克拉默法则相关。至于隐函数法则等等，不再本文的讨论范围内，有兴趣自行学习。
　　矩阵论的老师喜欢强调：行列式是没有意义的。虽然博主觉得老师想强调的是几何意义，但是行列式在几何的多用于做一些简单的判断，比如不同系统的变化趋势正负、线性相关、容积变化等。
　　行列式是解析数学中一个重要工具——克拉默法则解方程组。
　
　　设定两个可微的隐函数
$\begin{cases} F(x,y,z,u,v)=0\ G(x,y,z,u,v)=0 \end{cases}$
　　且 $u,v$ 可以用 $x,y,z$ 表示，即
$\begin{cases} F(x,y,z,f(x,y,z),g(x,y,z))=0\ G(x,y,z,f(x,y,z),g(x,y,z))=0 \end{cases}$
　　得到偏微分方程
$\begin{cases} \frac{\partial F}{\partial x}+\frac{\partial F}{\partial u}\frac{\partial u}{ \partial x}+\frac{\partial F}{\partial v}\frac{\partial v}{\partial x} &=&0 \ \frac{\partial G}{\partial x}+\frac{\partial G}{\partial u}\frac{\partial u}{ \partial x}+\frac{\partial G}{\partial v}\frac{\partial v}{\partial x} &=&0 \end{cases}$
　　使用克拉默法则解得
$\frac{\partial u}{\partial x}=-\frac{\det \begin{pmatrix} \frac{\partial F}{\partial x}& \frac{\partial F}{\partial v}\\ \frac{\partial G}{\partial x}& \frac{\partial G}{\partial v} \end{pmatrix}}{\det \begin{pmatrix} \frac{\partial F}{\partial u}& \frac{\partial F}{\partial v}\\\frac{\partial G}{\partial u}& \frac{\partial G}{\partial v} \end{pmatrix}}$
　　如果对行列式用一个符号来表示
$\frac{\partial (F,G)}{\partial (u,v)}=\det \begin{pmatrix} \frac{\partial F}{\partial u}&\frac{\partial F}{\partial v}\\ \frac{\partial G}{\partial u} &\frac{\partial G}{\partial v} \end{pmatrix}$
　　所以，可以解得
$\dfrac{\partial u}{\partial x}=-\dfrac{\dfrac{\partial (F,G)}{\partial (x,v)}}{\dfrac{\partial (F,G)}{\partial (u,v)}}$
　　这里 ${\dfrac{\partial (F,G)}{\partial (x,y)}},{\dfrac{\partial (F,G)}{\partial(u,v)}}$ 就被称为雅可比式（Jacobians)。

附：计算 $z$ 偏导数（梯度）的 GLSL 代码

vec2 DepthGradient(vec2 uv, float z)
{
  vec2 dz_duv = vec2(0.0, 0.0);

  vec3 duvz_dx = dFdx(vec3(uv,z));
  vec3 duvz_dy = dFdy(vec3(uv,z));

  dz_duv.x = duvz_dy.y * duvz_dx.z;
  dz_duv.x -= duvz_dx.y * duvz_dy.z;

  dz_duv.y = duvz_dx.x * duvz_dy.z;
  dz_duv.y -= duvz_dy.x * duvz_dx.z;

  float det = (duvz_dx.x * duvz_dy.y) - (duvz_dx.y * duvz_dy.x);
  dz_duv /= det;

  return dz_duv;
}

参考：
[1] Jacobian matrix and determinant – Wikipedia
[2] What is Jacobian Matrix?
[3] Cascaded Shadow Maps