Truncated Cauchy Non-negative Matrix Factorization

To analyze the generalization ability of Truncated CauchyNMF, we further assume that samples $[v_{1},\dots,v_{n}]$ are independent and identically distributed and drawn from a space $\mathcal{V}$ with a Borel measure $\rho$. We use $A_{\cdot j}$ and $A_{ij}$ to denote the $j$-th column and the $(i,j)$-th entry of a matrix , respectively, and $a_{i}$ is the $i$-th entry of a vector $a$.

For any $W\in\mathbb{R}_{+}^{m\times r}$, we define the reconstruction error of a sample $v$ as follows:

Therefore, the objective function of Truancated CauchyNMF (2) can be written as

Let us define the empirical reconstruction error of Truncated CauchyNMF as $R_{n}(f_{W})=\frac{1}{n}\sum_{i=1}^{n}f_{W}(v_{i})$, and the expected reconstruction error of Truncated CauchyNMF as $R(f_{W})=E_{v}\frac{1}{n}\sum_{i=1}^{n}f_{W}(v_{i})$. Intuitively, we want to learn

However, since the distribution of $v$ is unknown, we cannot minimize $R(f_{W})$ directly. Instead, we use the empirical risk minimization (ERM, [2]) algorithm to learn $W_{n}$ to approximate $W_{*}$, as follows:

We are interested in the difference between $W_{n}$ and $W_{*}$. If the distance is small, we can say that $W_{n}$ is a good approximation of $W_{*}$. Here, we measure the distance of their reduced expected reconstruction error as follows:

where $F_{\mathcal{W}}=\{f_{W}|W\in\mathcal{W}=\mathbb{R}_{+}^{m\times r}\}$. The right hand side is known as the generalization error. Note that since NMF is convex with respect to either $W$ or $H$ but not both, the minimizer $f_{W_{n}}$ is hard to obtain. In practice, a local minimizer is used as an approximation. Measuring the distance between the local minimizer and the global minimizer is also an interesting and challenging problem.

By analyzing the covering number [30] of the function class $F_{\mathcal{W}}$ and Lemma 1, we derive a generalization error bound for Truncated CauchyNMF as follows:

Remark 1. Theorem 1 shows that under the setting of our proposed Truncated CauchyNMF, the expected reconstruction error $R(f_{W_{n}})$ will converge to $R(f_{W_{*}})$ with a fast rate of order $O(\sqrt{{\ln n}/{n}})$, which means that when the sample size $n$ is large, the distance between $R(f_{W_{n}})$ and $R(f_{W_{*}})$ will be small. Moreover, if $n$ is large and a local minimizer $W$ (obtained by optimizing the non-convex objective of Truncated CauchyNMF) is close to the global minimizer $W_{n}$, the local minimizer will also be close to the optimal $W_{*}$.

Remark 2. Theorem 1 also implies that for any $W$ learned from (2), the corresponding empirical reconstruction error $R_{n}(f_{W})$ will converge to its expectation with a specific rate guarantee, which means our proposed Truncated CauchyNMF can generalize well to unseen data.

Note that the noise sampled from the Cauchy distribution should not be bounded because Cauchy distribution is heavy-tailed. And bounded observations always imply bounded noise. However, Theorem 1 keeps the boundedness assumption on the observations for two reasons: (1) the truncated loss function indicates that the observations corresponding to unbounded noise are discarded, and (2) in real applications, the energy of observations should be bounded, which means their $L_{2}$-norms are bounded.

We next compare the robustness of Truncated CauchyNMF with those of other NMF models by using a sample-weighted procedure interpretation [20]. The sample-weighted procedure compares the robustness of different algorithms from the optimization viewpoint.

Let $F(WH)$ denote the objective function of any NMF problem and $f(t)=F(tWH)$ where $t\in\mathbb{R}$. We can verify that the NMF problem is equivalent to finding a pair of $WH$ such that $f^{{}^{\prime}}(1)=0$³³3When minimizing $F(WH)$, the low rank matrices $W$ and $H$ will be updated alternately. Fixing one of them and optimizing the other implies that $f^{{}^{\prime}}(1)=0$. In other words, if $f^{{}^{\prime}}(1)\neq 0$, neither $W$ nor $H$ can be a minimizer., where $f^{{}^{\prime}}(t)$ denotes the derivative of $f(t)$. Let $c(V_{ij},WH)=(V-WH)_{ij}(-WH)_{ij}$ be the contribution of the $j$-th entry of the $i$-th training example to the optimization procedure and $e(V_{ij},WH)=|V-WH|_{ij}$ be an error function. Note that we choose $c(V_{ij},WH)$ as the basis of contribution because we choose NMF, which aims to find a pair of $WH$ such that $\sum_{ij}c(V_{ij},WH)=0$ and is sensitive to noise, as the baseline for comparing the robustness. Also note that $e(V_{ij},WH)$ represents the noise added to the $(i,j)$-th entry of $V$. The interpretation of the sample-weighted procedure explains the optimization procedure as being contribution-weighted with respect to the noise.

We compare $f^{\prime}(1)$ of a family of NMF models in Table I. Note that since multiplying $f^{\prime}(1)$ by a constant will not change its zero points, we can normalize the weights of different NMF models to unity when the noise is equal to zero. During the optimization procedures, robust algorithms should assign a small weight to an entry of the training set with large noise. Therefore, by comparing the derivative $f^{\prime}(1)$, we can easily make the following statements: (1) $L_{1}$-NMF⁴⁴4For the soundness of defining the subgradient of $L_{1}$-norm, we state that $\frac{0}{0}$ can be any value in $[-1,1]$. is more robust to noise and outliers than NMF; Huber-NMF combines the ideas of NMF and $L_{1}$-NMF; (2) HCNMF, $L_{2,1}$-NMF, RCNMF, and RNMF-$L_{1}$ work similarly to $L_{1}$-NMF because their weights are of order $O({1}/{e(V_{ij},WH)})$ with respect to the noise. It also becomes clear that HCNMF, $L_{2,1}$-NMF, and RCNMF exploit some data structure information because the weights include the neighborhood information of $e(V_{ij},WH)$ and that RNMF-$L_{1}$ is less sensitive to noise because it employs a sparse matrix $S$ to adjust the weights; (3) The interpretation of the sample-weighted procedure also illustrates why CIM-NMF works well for heavy noise. This is because its weights decrease exponentially when the noise is large; And (4) for the proposed Truncated CauchyNMF, when the noise is larger than a threshold, its weights will drop directly to zero, which decrease far faster than that of CIM-NMF and thus Truncated CauchyNMF is very robust to extreme outliers. Finally, we conclude that Truncated CauchyNMF is more robust than any other NMF models with respect to extreme outliers because it has the power to provide smaller weights to examples.

Generally speaking, the half-quadratic (HQ) programming algorithm [9] reformulates the non-quadratic loss function as an augmented loss function in an enlarged parameter space by introducing an additional auxiliary variable based on the convex conjugation theory [3]. HQ is equivalent to the quasi-Newton method [24] and has been widely applied in non-quadratic optimization.

Note that the function $f(x):\mathbb{R}_{+}\rightarrow\mathbb{R}$ is continuous, and according to [3], its conjugate $f^{*}(y):\mathbb{R}\rightarrow\mathbb{R}\cup\{+\infty\}$ is defined as

Since $f(x)$ is convex and closed (although the domain $\mathbb{R}_{+}$ is open, $f(x)$ is closed, see Section A.3.3 in [3]), the conjugate of its conjugate function is itself [3], i.e., $f^{**}=f$, then we have:

By substituting $x=(\frac{V-WH}{\gamma})_{ij}^{2}$ into (9), we have the augmented loss function

where $Y_{ij}$ is the auxiliary variable introduced by HQ for $(\frac{V-WH}{\gamma})_{ij}^{2}$. By substituting (10) into (8), we have the objective function in an enlarged parameter space

where the equality comes from the separability of the optimization problems with respect to $Y_{ij}$.

Although the objective function in (8) is non-quadratic, its equivalent problem (3.1) is essentially a quadratic optimization. In this paper, HQ solves (3.1) based on the block coordinate descent framework. In particular, HQ recursively optimizes the following three problems. At $t$-th iteration,

Using Theorem 2, we know that the solution of (12) can be expressed analytically as

Since (13) and (14) are symmetric and intrinsically weighted non-negative least squares (WNLS) problems, they can be optimized in the same way using the Nesterov method [10]. Taking (13) as an example, the procedure of its Nesterov based optimization is summarized in Algorithm 1, and its derivative is derived in the supplementary material. Considering that (13) is a constrained optimization problem, similar to [18], we use the following projected gradient-based criterion to check the stationarity of the search point , i.e., $\nabla_{j}^{P}(h_{k})=0$, where $\nabla_{j}^{P}(h_{k})_{l}=\left\{\begin{array}[]{cc}\nabla_{j}^{P}(h_{k})_{l},&(h_{k})_{l}\geq 0\\ \min\{0,\nabla_{j}^{P}(h_{k})_{l}\},&(h_{k})_{l}=0\end{array}\right.$. Since the above stopping criterion will make OGM run unnecessarily long, similar to [18], we use a relaxed version

where $\epsilon_{1}$ is a tolerance that controls how far the search point is from a stationary point.

The complete procedure of the HQ algorithm is summarized in Algorithm 2. The weights of entries and factor matrices are updated recursively until the objective function does not change. We use the following stopping criterion to check the convergence in Algorithm 2:

where $\epsilon_{2}$ signifies the tolerance, $F(W,H)$ signifies the objective function of (8) and $(W^{*},H^{*})$ signifies a local minimizer⁵⁵5Since any local minimal is unknown beforehand, we instead utilize $(W^{t-1},H^{t-1})$ in our experiments.. The stopping criterion (16) implies that HQ stops when the search point is sufficiently close to the minimizer and sufficiently far from the initial point. Line $3$ updates the scale parameter by the Nagy algorithm and will be further presented in Section 3.2. Line $4$ detects outliers by robust statistics and will be presented in Section 3.3.

The main time cost of Algorithm 2 is incurred on lines $2$, $4$, $5$, $6$, $7$, $8$, and $9$. The time complexities of lines $2$ and $7$ are both $O(mnr)$. According to Algorithm 1, the time complexities of lines $6$ and $9$ are $O(mr^{2})$ and $O(nr^{2})$, respectively. Since line $4$ introduces a median operator, its time complexity is $O(mn\ln(mn))$. In summary, the total complexity of Algorithm 2 is $O((mn\ln(mn)+mnr^{2}))$.

The parameter estimation problem for Cauchy distribution has been studied for several decades [32][33][34]. Nagy [32] proposed an I-divergence based method, termed the Nagy algorithm for short, to simultaneously estimate location and scale parameters. The Nagy algorithm minimizes the discrimination information⁶⁶6The discrimination information of random variable $\xi_{1}$ given random variable $\xi_{2}$ is defined as $D(\xi_{1}|\xi_{2})=\int_{-\infty}^{+\infty}\ln\frac{f_{1}(x)}{f_{2}(x)}dF_{1}(x)$, where $f_{1}$ and $f_{2}$ are the PDFs of $\xi_{1}$ and $\xi_{2}$, and $F_{1}$ is the distribution function of $\xi_{1}$. between the empirical distribution of the data points and the prior Cauchy distribution with respect to the parameters. In our Truncated CauchyNMF model (2), the location parameter of the Cauchy distribution is assumed to be zero, and thus we only need to estimate the scale-parameter $\gamma$.

Here we employ the Nagy algorithm to estimate the scale-parameter based on all the residual errors of the data. According to [32], supposing there exist a large number of residual errors, the scale-parameter estimation problem can be formulated as

where $D(\cdot|\cdot)$ denotes the discrimination information, and the first equality is due to the independence of $\eta_{n}$ and $\gamma$, and the second equality is due to the Law of large numbers. By substituting the probability density function $f_{0,\gamma}$ of Cauchy distribution⁷⁷7The probability density function (PDF) of Cauchy distribution is $f_{x_{0},\gamma}(x)={1}/({\pi\gamma(1+(\frac{x-x_{0}}{\gamma})^{2})})$, where $x_{0}$ is the location parameter, specifying the location of the peak of the distribution, and $\gamma$ is the scale parameter, specifying the half-width at half-maximum. into (17) and replacing $\{x_{n}\}$ with $\{E_{ij}\}$, we can rewrite (17) as follows: $\min_{\gamma}\sum_{i=1}^{m}\sum_{j=1}^{n}\frac{1}{mn}\ln\{\pi\gamma(1+(\frac{E_{ij}}{\gamma})^{2})\}$. To solve this problem, Nagy [32] proposed an efficient iterative algorithm, i.e.,

where $\gamma_{0}>0$, and $e_{k}^{0}=\frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}\frac{1}{(1+(\frac{E_{ij}}{\gamma_{k}})^{2})}$. In [32], Nagy proved that the algorithm (18) converges to a fixed point assuming the number of data points is large enough, and this assumption is reasonable in Truncated CauchyNMF.

Looking more carefully at (12), (13) and (14), HQ intrinsically assigns a weight for each entry of $V$ with both factor matrices $H^{t+1}$ and $W^{t+1}$ fixed, i.e., $Q_{ij}^{t+1}=\left\{\begin{array}[]{cc}\frac{1}{1+(\frac{E^{t}}{\gamma})_{ij}^{2}},&if\;|E^{t}_{ij}|\leq\gamma\sqrt{\sigma}\\ 0,&if\;|E^{t}_{ij}|>\gamma\sqrt{\sigma}\end{array}\right.$, where $E^{t}$ denotes the error matrix at the $t$-th iteration. The larger the magnitude of error for a particular entry, the lighter the weight is assigned to it by HQ. Intuitively, the corrupted entry contributes less in learning the intrinsic subspace. If the magnitude of error exceeds a threshold $\gamma\sqrt{\sigma}$, Truncated CauchyNMF assigns zero weights to the corrupted entries to inhibit their contribution to the learned subspace. That is how Truncated CauchyNMF filters out extreme outliers.

However, it is non-trivial to estimate the threshold $\gamma\sqrt{\sigma}$. Here, we introduce a robust statistics-based method to explicitly detect the support of the outliers instead of estimating the threshold to detect outliers. Since the energy function of Truncated CauchyNMF gets close to that of NMF as the error tends towards zero, i.e., $\lim_{x\rightarrow 0}(\ln(1+x^{2})-x^{2})=0$. Truncated CauchyNMF encourages the small magnitude errors to have a Gaussian distribution. Let $\Theta^{t}$ denote the set of magnitudes of error at the $t$-th iteration of HQ, i.e., $\Theta^{t}=\{|E_{ij}^{t}|:1\leq i\leq m,1\leq j\leq n\}$ where $E^{t}=V-W^{t}H^{t}$. It is reasonable to believe that a subset of $\Theta^{t}$, i.e., $\Gamma^{t}=\{\theta\in\Theta^{t}:\theta\leq med\{\Theta^{t}\}\}$, obeys a Gaussian distribution, where $med\{\Theta^{t}\}$ signifies the median of $\Theta^{t}$. Since $|\Gamma^{t}|=\lfloor\frac{mn}{2}\rfloor$, it suffices to estimate both the mean $\mu^{t}$ and standard deviation $\delta^{t}$ from $\Gamma^{t}$. According to the three-sigma-rule, we detect the outliers as $O^{t}=\{\tau\in\Gamma^{t}:|\tau-\mu^{t}|>3\delta^{t}\}$ and output their indices $\Omega(t)$.

To illustrate the effect of outlier rejection, Figure 3 presents a sequence of weighting matrices generated by HQ for the motivating example described in Figure 2. It shows that HQ correctly assigns zero weights for the corrupted entries in only a few iterations and finally detects almost all outliers including illumination, sunglasses, and scarves (see the last column in Figure 3) in the end.

Traditional NMF [17] assumes that noise obeys a Gaussian distribution and derives the following squared $L_{2}$-norm based objective function: $\min_{W\geq 0,H\geq 0}\|V-WH\|_{F}^{2}$, where $\|X\|_{F}=\sqrt{\sum_{ij}X_{ij}^{2}}$ signifies the matrix Frobenius norm. It is commonly known that NMF can be solved by using the multiplicative update rule (MUR, [17]). Because of the nice mathematical property of squared $L_{2}$-norm and the efficiency of MUR, NMF has been extended for various applications [4][6][28]. However, NMF and its extensions are non-robust because the $L_{2}$-norm is sensitive to outliers.

Hamza and Brady [12] proposed a hypersurface cost based NMF (HCNMF) by minimizing the summation of hypersurface costs of errors, i.e., $\min_{W\geq 0,H\geq 0}\{\sum_{ij}\delta((V-WH)_{ij})\}$, where $\delta(x)=\sqrt{1+x^{2}}-1$ is the hypersurface cost function. According to [12], the hypersurface cost function has differentiable and bounded influence function. Since the hypersurface cost function is differentiable, HCNMF can be directly solved by using the projected gradient method. However, the optimization of HCNMF is difficult because Armijo’s rule based line search is time consuming [12].