## Nonnegative Matrix Factorization


```python
%matplotlib inline
%load_ext Cython

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import NMF as NMF_sklearn
```

NMF is a method to factorize a nonnegative matrix $X_{m \times n}$ to nonnegative matrix factors $W_{m \times r}$ and $H_{r \times n}$ such that $X \approx WH$. Here, $r$ is a parameter to be determined and can be called the inner dimension.

### Iterative Algorithm for NMF
An iterative algorithm for NMF as proposed in [[1]](#references) using the multiplicate update rules. This implementation is only for practice; when you need to use NMF, please refer to [sklearn.decomposition.NMF](http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.NMF.html).


```cython
%%cython
import numpy as np
cimport numpy as np
cimport cython

DTYPE = np.float64
ctypedef np.float_t DTYPE_t

@cython.boundscheck(False)
@cython.wraparound(False)
def obj_func(np.ndarray[DTYPE_t, ndim=2] orig, np.ndarray[DTYPE_t, ndim=2] approx):
    """
    Compute the objective function value for NMF given the original
    and approximate matrices.
    
    Parameters
    ----------
    orig : numpy.array
        Original matrix to factorize.
    approx : numpy.array
        Current approximate matrix.
        
    Returns
    -------
    obj_value : float
        Current objective function value.
    """
    return np.sum(orig*np.log(approx) - approx)

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def nmf(np.ndarray[DTYPE_t, ndim=2] X, int r):
    """
    Given a nonnegative matrix X of shape m x n, iteratively
    compute a nonnegative matrix factorization X = WH, and
    return the factors W and H.
    
    The returned factors have shapes (m, r) and (r, n),
    respectively.
    
    Parameters
    ----------
    X : numpy.array
        m x n nonnegative matrix.
    r : int
        Number of columns of factor W.
        
    Returns
    -------
    W : numpy.array
        m x r nonnegative matrix.
    H : numpy.array
        r x n nonnegative matrix.
    """
    cdef int m = X.shape[0]
    cdef int n = X.shape[1]
    cdef int prev_index = 0
    cdef int a, i, mu
    cdef DTYPE_t curr_obj, curr_mean
    cdef DTYPE_t prev_mean = 0
    cdef DTYPE_t coeff, col_sum
    cdef np.ndarray[DTYPE_t, ndim=2] W = np.random.rand(m, r)
    cdef np.ndarray[DTYPE_t, ndim=2] H = np.random.rand(r, n)
    cdef np.ndarray[DTYPE_t, ndim=1] prev_objs = np.zeros(10, dtype=DTYPE)
    cdef np.ndarray[DTYPE_t, ndim=2] XA
    
    while True:
        XA = W@H 
        
        # Convergence check
        curr_obj = obj_func(X, XA)
        if abs(curr_obj - prev_mean) < 1e-6:
            break
        else:
            prev_mean = (10*prev_mean - prev_objs[prev_index] + curr_obj)/10
            prev_objs[prev_index] = curr_obj
            prev_index = (prev_index + 1)%10
            
        # Update each element of W while fixing H
        for i in range(m):
            for a in range(r):
                coeff = 0
                for mu in range(n):
                    coeff += X[i, mu]/XA[i, mu]*H[a, mu]
                W[i, a] *= coeff
        
        # Normalize columns of W st they sum to 1
        for a in range(r):
            col_sum = 0
            for i in range(m):
                col_sum += W[i, a]
            for i in range(m):
                W[i, a] /= col_sum
        
        # Update each element of H while fixing W
        for a in range(r):
            for mu in range(n):
                coeff = 0
                for i in range(m):
                    coeff += W[i, a]*X[i, mu]/XA[i, mu]
                H[a, mu] *= coeff
                
    return W, H
```


```python
X = np.random.randint(low=0, high=100, size=(10, 200)).astype(np.float)
```


```python
r = 4
W, H = nmf(X, r)

learner = NMF_sklearn(n_components=r, solver='mu', max_iter=1000)
W_sklearn = learner.fit_transform(X)
H_sklearn = learner.components_
```


```python
D = X - W@H
D_sklearn = X - W_sklearn@H_sklearn

print('Divergence of our impl.    : {}'.format(la.norm(D, ord='fro')))
print('Divergence of sklearn impl.: {}'.format(la.norm(D_sklearn, ord='fro')))
```

    Divergence of our impl.    : 976.5349142405947
    Divergence of sklearn impl.: 1009.184869947594


## Bayesian Inference For NMF

### $r = 1$ Case
In this case, we factorize $X$ into row and column vectors $W$ and $H$ whose outer product is the approximation $\tilde{X}$.

$$
\begin{align}
    \begin{pmatrix}
      x_{11} & \dots & x_{1n} \\
      \vdots & \ddots & \vdots \\
      x_{m1} & \dots & x_{mn}
    \end{pmatrix}
    \approx
    \begin{pmatrix}
      w_1 \\
      \vdots \\
      w_m
    \end{pmatrix}
    \begin{pmatrix}
      h_1 & \dots & h_n
    \end{pmatrix}
\end{align}
$$


```python
def rank_one_nmf(X):
    """
    Given a nonnegative matrix X, construct nonnegative vectors
    w and h whose outer product approximates X.
    """
    m, n = X.shape
    w = X[:, 0]
    h = np.ones(n)
    for i in range(1, n):
        w_unit = w/la.norm(w)
        proj = np.dot(X[:, i], w_unit)*w_unit   
        h[i] = np.mean(proj/w)
    return w, h
```


```python
w = np.array([5, 3, 7])
h = np.array([0.1, 2, 3, 0.7])
X = np.outer(w, h)
#X = np.random.randint(low=0, high=100, size=(3, 3))
print('Original matrix\n', X)
w_found, h_found = rank_one_nmf(X)
print('Reconstructed matrix\n', np.outer(w_found, h_found))
print('w = ', w_found)
print('h = ', h_found)
```

    Original matrix
     [[  0.5  10.   15.    3.5]
     [  0.3   6.    9.    2.1]
     [  0.7  14.   21.    4.9]]
    Reconstructed matrix
     [[  0.5  10.   15.    3.5]
     [  0.3   6.    9.    2.1]
     [  0.7  14.   21.    4.9]]
    w =  [ 0.5  0.3  0.7]
    h =  [  1.  20.  30.   7.]


<a id='references'></a>
## References
1. D. Lee and H. S. Seung, “Learning the parts of objects with nonnegative matrix factorization,” Nature, vol. 401, pp. 788– 791, 1999.

e9a71485cfd3f67c3a46f076bb449476a7f4dac1
ipynb
Jupyter Notebook
nonnegative_matrix_factorization.ipynb
eozd/experiments
8974235c9cf6610e7c576c04b4cd50f4f421b2a7
MIT
nonnegative_matrix_factorization.ipynb
eozd/experiments
8974235c9cf6610e7c576c04b4cd50f4f421b2a7
MIT
nonnegative_matrix_factorization.ipynb
eozd/experiments
8974235c9cf6610e7c576c04b4cd50f4f421b2a7
MIT
Qwen/Qwen-72B
1. YES
2. YES