在任意向量的正交空间中采样

有时候我们需要在特定的空间中采样，观察样本的性质。

现有一组$D$维空间中的向量${\mathbf{x}}^1,{\mathbf{x}}^2,\cdots {\mathbf{x}}^N$，其中${\mathbf{x}}^n=(x_1^n,x_2^n,\cdots ,x_D^n)$，($D>N$)，采样$\mathbf{s}$，使得$\mathbf{s}$与所有${\mathbf{x}}^n$正交。

我们要找$x^n$空间的一组基底，用特征值分解就好。$$X=U\Sigma V^{\top }$$
其中，$X$是按行排列的 x n , n ∈ { 1 , 2 , ⋯ N } \mathbf{x}^n, n\in\{1,2,\cdots N \} xn,n∈{1,2,⋯N}，并用$0$填充空行。
$$\Sigma =\left(\begin{array}{cccccccc}{\lambda }_1& & & & & & & \\ & {\lambda }_2& & & & & & \\ & & \ddots & & & & & \\ & & & {\lambda }_K& & & & \\ & & & & 0& & & \\ & & & & & \ddots & & \\ & & & & & & 0& \end{array}\right)$$
$K$ 为非零特征值个数。
即通过线性变换$U$和$V^{\top }$，我们得到${\mathbf{x}}^n$的一组基底$\Sigma$。显然，与$\Sigma$正交的空间的一组基底为
$$B=\left(\begin{array}{cc}{O}_{K\times K}& O_{K\times (D-K)}\\ O_{(D-K)\times K}& I_{(D-K)\times (D-K)}\end{array}\right)$$
$O$为零矩阵，$I$为单位矩阵。那么，与${\mathbf{x}}^n$正交的空间的一组基底为$UB{V}^{\top }$。

import numpy as np

z_dim = 512 # number of dimensions
n_vectors = 10 # number of Xs
assert n_vectors < z_dim
n_samples = 1000 # number of samples to orthogonal Xs

Xs = np.random.rand(n_vectors, z_dim)
Xs = np.concatenate((Xs,np.zeros((z_dim-n_vectors, z_dim))))

U, s, V = np.linalg.svd(Xs)
assert np.linalg.norm(np.dot(np.dot(U, np.diag(s)), V) - Xs) < 1e-8

num_lambdas = np.sum(s>1e-8)

bases = np.eye(z_dim)
bases[range(num_lambdas),range(num_lambdas)] = 0
bases = np.dot(np.dot(U, bases), V)

samples = np.random.rand(z_dim, n_samples)
samples = np.dot(bases.T, samples)

assert np.mean(np.abs(np.dot(Xs, samples))) < 1e-8

另外，还可以使用旋转矩阵，参见python n维向量向任意方向旋转，求旋转矩阵

= = =
正文到此为止，下面是一些小的尝试。

接下来，给定目标向量$x^t$，可以计算其在正交空间中的投影，正交化，或者，只靠近$x^t$而与 x i , i ∈ { 1 , 2 , ⋯ N } x^i, i\in \{1,2,\cdots N \} xi,i∈{1,2,⋯N}距离不变。

x_target = np.random.rand(1, z_dim)
x_proj = np.concatenate((x_target, np.zeros((z_dim-1,z_dim))))
x_proj = np.dot(np.dot(U.T, x_proj), V.T)
# x_proj = x_proj[:,0]
x_proj[:, :num_lambdas] = 0
x_proj = np.dot(np.dot(U, x_proj), V)

assert np.linalg.norm(np.dot(Xs, x_proj[0,:])) < 1e-8

# The defination of projection matrix
P = np.dot(Xs.T, np.linalg.inv((Xs + bases).T))
x_proj_check = np.dot(np.eye(z_dim)-P, x_target[0,:])

assert np.linalg.norm(np.dot(Xs, x_proj_check)) < 1e-8
assert np.linalg.norm(x_proj[0,:]-x_proj_check) < 1e-8

从原点出发，向$x^t$移动而与 x i , i ∈ { 1 , 2 , ⋯ N } x^i, i\in \{1,2,\cdots N \} xi,i∈{1,2,⋯N}相切的点的运动轨迹

import numpy as np

z_dim = 512 # number of dimensions
n_vectors = 2 # number of Xs
assert n_vectors < z_dim

Xs = np.random.rand(n_vectors, z_dim)
Xs = np.concatenate((Xs,np.zeros((z_dim-n_vectors, z_dim))))

x_target = np.random.rand(1, z_dim)

# debug for n_vector=1
dist_X_target = np.linalg.norm(Xs[:n_vectors,:]-x_target, axis=1)
# debug for n_vector=1

while True:
        U, s, V = np.linalg.svd(Xs)

        num_lambdas = np.sum(s>1e-8)

        bases = np.eye(z_dim)
        bases[range(num_lambdas),range(num_lambdas)] = 0
        bases = np.dot(np.dot(U, bases), V)

        x_proj = np.concatenate((x_target, np.zeros((z_dim-1,z_dim))))
        x_proj = np.dot(np.dot(U.T, x_proj), V.T)
        # x_proj = x_proj[:,0]
        x_proj[:, :num_lambdas] = 0
        x_proj = np.dot(np.dot(U, x_proj), V)

        assert np.linalg.norm( np.dot(Xs, x_proj[0,:])) < 1e-8

        x_step = 0.1 * x_proj[[0],:]
        Xs[:n_vectors,:] = Xs[:n_vectors,:] - x_step
        x_target = x_target - x_step

        dist_Xs_step = np.linalg.norm(Xs)
        dist_target_step = np.linalg.norm(x_target)
        # print([dist_Xs_step, dist_target_step])
        dist = np.linalg.norm(x_target),np.linalg.norm(Xs[:n_vectors,:], axis=1)
        print([dist_X_target, dist])
        if np.linalg.norm(x_step)<1e-5:
                break

上面假设了$D>N$。当$D\le N$，且$rank(X)<D$，我们仍可以用上面的方法。当$X$张成整个空间，即$rank(X)=D$，无法找到一个切线，这是一个只有零解的线性规划问题：
$$\begin{gathered} \underset{v}{\max }{x}^{t}v \\ s.t.Xv=0 \end{gathered}$$
如果将这个问题进行一下放缩，例如用拉格朗日乘子法：$$\underset{v,\delta }{\max }{x}^{t}v-\delta |Xv|$$
$\delta \in {\mathbb{R}}^N$是拉格朗日乘子，还可以继续求解。这里不展开了。

当然，还有一个取巧的办法，仍然使用svd分解。矩阵$\Sigma$中，特征值按从大到小的顺序排列。特征值的大小其实代表了这一维度的权重，最大的特征值对应的主方向，代表样本朝向这个主方向的分量最多（参考主成分分析PCA）。那么我们取特征值最小的主方向，沿着这个方向投影$x^t$即得所求。