matching pursuit笔记

给定一个 Hilbert space 中的向量 $f$ 和一个称为字典的向量组 D = { x 1 , … , x n } D=\{x_1,\dots,x_n\} D={x1​,…,xn​}，且 $\parallel x_i\parallel =1$，想用 $D$ 中的某些元素的组合尽量逼近 $f$，反过来说，是将 $f$ 拆成 $D$ 中元素的组合。

Matching Pursuit

MP [1] 的思路是递归地增量逼近。当随便拿一个 $x_i$ 近似时，MP 只取 $f$ 在 $x_i$ 上正交投影，此时 f 可以拆成：$$\begin{array}{rl}f& =\hat{f}+Rf\\ & =<f,x_i>x_i+Rf\end{array}$$ 其中 $\hat{f}$ 是近似值，$<\cdot ,\cdot >$ 是内积，$Rf$ 表示 $\hat{f}$ 与 $f$ 之间的误差。

为了使近似效果最好，应选一个使得 $|<f,x_i>|$ 最大的 $x_i$，[1] 中说有时需要将这个限制松弛成 $$|<f,x_i>|\ge \alpha \cdot \underset{j}{\sup }|<f,x_j>|,0<\alpha <1$$

选定上述的最优 $x_i$ 之后，可以继续用 $D$ 拆误差 $Rf$：$$Rf=<Rf,x_{i^{\prime }}>x_{i^{\prime }}+R(Rf)$$ 于是原来的近似可以修正成：$$\begin{array}{rl}f& =<f,x_i>x_i+Rf\\ & =\underset{\hat{f}}{\underbrace{<f,x_i>x_i+<Rf,x_{i^{\prime }}>x_{i^{\prime }}}}+R(Rf)\end{array}$$ 此时有 $$\parallel Rf{\parallel }^2=|<Rf,x_{i^{\prime }}>{|}^2+\parallel R(Rf){\parallel }^2\ge \parallel R(Rf){\parallel }^2$$ 即近似误差更小了。这可以导出一个迭代算法，表示成：$$f\approx {\hat{f}}_m=\sum _{i=1}^m<R_{i-1}f,x_{{\sigma }_i}>x_{{\sigma }_i}+R_{m}f$$ 其中 $x_{{\sigma }_i}$ 是第 i 步选的向量，$R_{i}f$ 是第 i 步的误差，$R_{0}f=f$。

记 V = span { x 1 , … , x n } ‾ V=\overline{\text{span}\{x_1,\dots,x_n\}} V=span{x1​,…,xn​}​ 是 $D$ 张成的空间（closed linear span of the dictionary vectors），则当 $V=H$ 时，$f$ 必能被 $D$ 组合得出即 $$\begin{gathered} \underset{m\to +\infty }{\lim }{\hat{f}}_m=f \\ \underset{m\to +\infty }{\lim }{R}_{m}f=0 \end{gathered}$$

注意，MP 中每个 $x_i$ 允许多次被选中。

Orthogonal Matching Pursuit

OMP [2] 表示 MP 只保证 $R_{i}f\perp x_{{\sigma }_i}$，想收敛可能需要很多步迭代。当限制迭代次数 $m$ 时，MP 得出的结果不是最优的。

OMP 增加了对 $R_{i}f$ 的约束，使其与所有 $x_{{\sigma }_j}(j=1,\dots ,i-1)$ 都正交（full backward orthogonality）。

假设在第 k 步的拆分为：$$f=\sum _{i=1}^k{a}_i^k{x}_{{\sigma }_i}+R_{k}f\text{\hspace{1em}(1)}$$ 其中 $a_i^k$ 是第 k 步中的组合系数，$R_{k}f$ 满足 $R_{k}f\perp x_{{\sigma }_j},j=1,\dots ,k$，即 full backward orthogonality。在第 k + 1 步，按 MP 的方式选 $x_{{\sigma }_{k+1}}$，新的拆分为：$$f=\sum _{i=1}^{k+1}{a}_i^{k+1}{x}_{{\sigma }_i}+R_{k+1}f\text{\hspace{1em}(2)}$$ 约束 $R_{k+1}f$ 同样满足 full backward orthogonality。为了确定系数 $a_i^{k+1}$，先将 $x_{{\sigma }_{k+1}}$ 同样用 { x σ 1 , … , x σ k } \{x_{\sigma_1},\dots,x_{\sigma_k}\} {xσ1​​,…,xσk​​} 拆成：$$x_{{\sigma }_{k+1}}=\sum _{i=1}^k{b}_i^k{x}_{{\sigma }_i}+{\gamma }_k\text{\hspace{1em}(3)}$$ 其中 $b_i^k$ 和 ${\gamma }_k$ 都是用某种方法算出来的。${\gamma }_k$ 是 $x_{{\sigma }_{k+1}}$ 中用 { x σ 1 , … , x σ k } \{x_{\sigma_1},\dots,x_{\sigma_k}\} {xσ1​​,…,xσk​​} 组合不出的成分，所以 ${\gamma }_k\perp x_{{\sigma }_j},j=1,\dots ,k$。将 (3) 代入 (2) 与 (1) 对比系数：$$\begin{array}{rl}f& =\underset{1}{\underbrace{\sum _{i=1}^k{a}_i^k{x}_{{\sigma }_i}}}+\underset{2}{\underbrace{R_{k}f}}\\ & =\sum _{i=1}^k{a}_i^{k+1}{x}_{{\sigma }_i}+a_{k+1}^{k+1}{x}_{{\sigma }_{k+1}}+R_{k+1}f\\ & =\sum _{i=1}^k{a}_i^{k+1}{x}_{{\sigma }_i}+a_{k+1}^{k+1}\left[\sum _{j=1}^k{b}_j^k{x}_{{\sigma }_j}+{\gamma }_k\right]+R_{k+1}f\\ & =\underset{1}{\underbrace{\sum _{i=1}^k[a_i^{k+1}+a_{k+1}^{k+1}\cdot b_i^k]x_{{\sigma }_i}}}+\underset{2}{\underbrace{a_{k+1}^{k+1}{\gamma }_k+R_{k+1}f}}\end{array}$$ 所以 $$\begin{gathered} a_i^{k+1}=a_i^k-a_{k+1}^{k+1}\cdot b_i^k \\ {R}_{k+1}f=R_{k}f-a_{k+1}^{k+1}\cdot {\gamma }_k \end{gathered}$$ 只要确定 $a_{k+1}^{k+1}$ 就可以。又由 $R_{k+1}f\perp x_{{\sigma }_{k+1}}$ 的限制，有：$$\begin{array}{rl}0& =<R_{k+1}f,x_{{\sigma }_{k+1}}>\\ & =<R_{k}f,x_{{\sigma }_{k+1}}>-a_{k+1}^{k+1}\cdot <{\gamma }_k,x_{{\sigma }_{k+1}}>\end{array}$$ 所以 $a_{k+1}^{k+1}=\frac{<R_{k}f,x_{{\sigma }_{k+1}}>}{<{\gamma }_k,x_{{\sigma }_{k+1}}>}$。

这样在保证 $R_{k+1}f\perp x_{{\sigma }_{k+1}}$ 的限制的同时，因为 $R_{k}f\perp x_{{\sigma }_j}\wedge {\gamma }_k\perp x_{{\sigma }_j}$，所以 $R_{k+1}f=R_{k}f-{\gamma }_k\perp x_{{\sigma }_j},j=1,\dots ,k$，所以 $R_{k+1}f$ 也满足 full backward orthogonality。

Non-negative Orthogonal Matching Pursuit

NOMP [3] 限制组合系数非负，说是可以提高稳定性。

[4] 是一份参考实现。下面是我改写的版本，跟 [4] 对拍过，结果一致。

# nomp.py
import numpy as np
import scipy.io as sio
import scipy.linalg
import scipy.optimize
from sklearn.preprocessing import normalize


"""Non-negative Orthogonal Matching Pursuit

ref:
- Stable and efficient representation learning with nonnegativity constraints
- Supplementary Materials: Compositional Zero-Shot Learning via Fine-Grained Dense Feature Composition
- https://github.com/hbdat/neurIPS20_CompositionZSL/blob/main/omp/omp.py#L43
"""


def to_pos(X):
    """X -> [max(0, X), max(0, -X)]
    把含有负值的向量映射成全是非负的向量
    相当于用多一倍的维度代偿原本负值的空间
    """
    X_inv = - X
    X_pos = np.concatenate([
        np.where(X > 0, X, 0),
        np.where(X_inv > 0, X_inv, 0),
    ], axis=1)
    print(X_pos.shape)
    return X_pos


def cos(A, B=None):
    """cosine"""
    An = normalize(A, norm='l2', axis=1)
    if (B is None) or (B is A):
        return np.dot(An, An.T)
    Bn = normalize(B, norm='l2', axis=1)
    return np.dot(An, Bn.T)


def spherical_kmeans(X, k=-1):
    """select a subset for dictionary using spherical k-means
    X: [n, d]
    k: # of cluster

    - https://github.com/jasonlaska/spherecluster
    """
    from spherecluster import SphericalKMeans
    n = X.shape[0]
    if (k <= 0) or (k > n):
        k = max(1, n // 2)
    skm = SphericalKMeans(n_clusters=k)
    skm.fit(X)
    return skm.cluster_centers_


def norm2(x):
    """x: [d] vector"""
    return np.linalg.norm(x) / np.sqrt(len(x))


def NOMP(D, f, k=-1, non_neg=True, max_iter=200, tol_res=1e-3, tol_coef=1e-12, verbose=False):
    """Non-negative Orthogonal Matching Pursuit
    Input:
        D: [m, d], dictionary (assumed to have been l2-normalised)
        f: [d], feature to decompose
        k: max # of selected component from dictionary,
            default to the whole dictionary size
        non_neg: to enforce non-negative coefficients
        tol_res: residual tolerance,
            exit if ||Rf|| < tol_res * ||f|| (error is small enough)
        tol_coef: residual covariance threshold,
            exit if coef_i < tol_coef * ||f|| for all i
    Output:
        active: [n] (n <= k), id of the selected components in the dictionary
        coef: [n], composition weights of the selected components
        Rf: [d], redidual vector
    """
    assert (2 == D.ndim) and (1 == f.ndim)
    assert D.shape[1] == f.shape[0], "dimension mismatch"
    m = D.shape[0]
    if (k <= 0) or (k > m):
        k = m

    f_norm = norm2(f)
    if verbose:
        print("initial residual:", f_norm)
    # tol_res = tol_res * f_norm
    tol_coef = tol_coef * f_norm
    Rf = f  # (initial) residual vector
    D_col = D.T  # [d, m], arrange as column vector
    active = []  # meta id of the selected components
    coef = np.zeros([m], dtype=np.float32)

    # degenerate solution
    if f_norm < tol_res:
        return active, coef

    if non_neg:
        get_candidate = lambda score: np.argmax(score)
        calc_coef = lambda _D, _f: scipy.optimize.nnls(_D, _f)[0]
    else:
        get_candidate = lambda score: np.argmax(np.abs(score))
        calc_coef = lambda _D, _f: scipy.linalg.lstsq(_D, _f)[0]

    for _it in range(max_iter):
        # print("---", _it, "---")
        # neglect the numerator cuz `Dict` is already normalised
        score = cos(Rf[np.newaxis, :], D)[0]  # [m]
        d_id = get_candidate(score)
        max_socre = score[d_id]
        if max_socre < tol_coef:
            break

        if d_id not in active:
            active.append(d_id)

        _coef = calc_coef(D_col[:, active], f)
        coef[active] = _coef

        f_hat = D_col[:, active].dot(_coef).flatten()
        Rf = f - f_hat

        err_rate = norm2(Rf) / f_norm
        # print("residual error rate:", err_rate)
        if err_rate < tol_res:
            break
        if len(active) >= k:
            break

    f_hat = D_col[:, active].dot(coef[active]).flatten()
    Rf = f - f_hat
    err = norm2(Rf)
    err_rate = err / f_norm
    if verbose:
    	# print("residual:", np.linalg.norm(res))
        print("final residule: {:.4f}, error rate: {:.4f}".format(err, err_rate))
    return active, coef, f_hat, Rf

References

1. Matching Pursuits With Time-Frequency Dictionaries
2. Orthogonal Matching Pursuit: Recursive Function Approximat ion with Applications to Wavelet Decomposition
3. Stable and Efficient Representation Learning with Nonnegativity Constraints
4. hbdat/neurIPS20_CompositionZSL/omp/omp.py
5. jasonlaska/spherecluster