Python Неотрицательная матричная факторизация, которая обрабатывает как нули, так и отсутствующие данные?

Используя этот лист преобразования кода Matlab в Python, я смог переписать NMF из Matlab toolbox.
Мне пришлось разложить матрицу 40k X 1k с разреженностью 0,7%. Используя 500 скрытых функций, моей машине потребовалось 20 минут на 100 итераций.

Вот метод:

import numpy as np
from scipy import linalg
from numpy import dot

def nmf(X, latent_features, max_iter=100, error_limit=1e-6, fit_error_limit=1e-6):
    """
    Decompose X to A*Y
    """
    eps = 1e-5
    print 'Starting NMF decomposition with {} latent features and {} iterations.'.format(latent_features, max_iter)
    X = X.toarray()  # I am passing in a scipy sparse matrix

    # mask
    mask = np.sign(X)

    # initial matrices. A is random [0,1] and Y is A\X.
    rows, columns = X.shape
    A = np.random.rand(rows, latent_features)
    A = np.maximum(A, eps)

    Y = linalg.lstsq(A, X)[0]
    Y = np.maximum(Y, eps)

    masked_X = mask * X
    X_est_prev = dot(A, Y)
    for i in range(1, max_iter + 1):
        # ===== updates =====
        # Matlab: A=A.*(((W.*X)*Y')./((W.*(A*Y))*Y'));
        top = dot(masked_X, Y.T)
        bottom = (dot((mask * dot(A, Y)), Y.T)) + eps
        A *= top / bottom

        A = np.maximum(A, eps)
        # print 'A',  np.round(A, 2)

        # Matlab: Y=Y.*((A'*(W.*X))./(A'*(W.*(A*Y))));
        top = dot(A.T, masked_X)
        bottom = dot(A.T, mask * dot(A, Y)) + eps
        Y *= top / bottom
        Y = np.maximum(Y, eps)
        # print 'Y', np.round(Y, 2)


        # ==== evaluation ====
        if i % 5 == 0 or i == 1 or i == max_iter:
            print 'Iteration {}:'.format(i),
            X_est = dot(A, Y)
            err = mask * (X_est_prev - X_est)
            fit_residual = np.sqrt(np.sum(err ** 2))
            X_est_prev = X_est

            curRes = linalg.norm(mask * (X - X_est), ord='fro')
            print 'fit residual', np.round(fit_residual, 4),
            print 'total residual', np.round(curRes, 4)
            if curRes < error_limit or fit_residual < fit_error_limit:
                break

return A, Y

Здесь я использовал разреженную матрицу Scipy в качестве входных данных, а пропущенные значения были преобразованы в 0 с использованием метода toarray(). Поэтому маска была создана с помощью функции numpy.sign(). Однако, если у вас есть значения nan, вы можете получить те же результаты, используя функцию numpy.isnan().

В Scipy есть метод для решить неотрицательную задачу наименьших квадратов (NNLS). В этом ответе я воспроизвожу свой сообщение в блоге об использовании NNLS scipy для не -отрицательная матричная факторизация. Вас также могут заинтересовать другие мои сообщения в блоге, в которых используется autograd, Tensorflow и CVXPY для NNMF.

Цель: нашей цели дана матрица A, разложите ее на два неотрицательных множителя следующим образом:

A (M×N) ≈ W (M×K) × H (K×N), such that  W (M×K) ≥ 0 and  H (K×N) ≥ 0 

Обзор

Наше решение состоит из двух шагов. Сначала мы фиксируем W и изучаем H при заданном A. Затем фиксируем H и изучаем W при заданном A. Мы повторяем эту процедуру итеративно. Фиксация одной переменной и изучение другой (в этом параметре) широко известна как чередующийся метод наименьших квадратов, поскольку проблема сводится к задаче наименьших квадратов. Однако важно отметить, что, поскольку мы хотим, чтобы ограничения W и H были неотрицательными, мы используем NNLS вместо метода наименьших квадратов.

Шаг 1: Изучение H, учитывая A и W

Используя приведенную выше иллюстрацию, мы можем узнать каждый столбец H, используя соответствующий столбец из A и матрицу W.

H[:,j]=NNLS(W,A[:,j])

Обработка отсутствующих записей в A

В задаче совместной фильтрации A обычно представляет собой матрицу пользовательских элементов, и в ней много пропущенных элементов. Эти отсутствующие записи соответствуют пользователям, которые не оценили элементы. Мы можем изменить нашу формулировку, чтобы учесть эти недостающие элементы. Учтите, что M' ≤ M записей в A имеют наблюдаемые данные, теперь мы должны изменить приведенное выше уравнение следующим образом:

H[:,j]=NNLS(W[mask],A[:,j][mask])

где маска находится путем рассмотрения только записей M '.

Шаг 2: Изучение W, учитывая A и H

Пример кода

Импорт

import numpy as np
import pandas as pd

Создание матрицы для факторизации

M, N = 20, 10

np.random.seed(0)
A_orig = np.abs(np.random.uniform(low=0.0, high=1.0, size=(M,N)))
print pd.DataFrame(A_orig).head()


0   1   2   3   4   5   6   7   8   9
0   0.548814    0.715189    0.602763    0.544883    0.423655    0.645894    0.437587    0.891773    0.963663    0.383442
1   0.791725    0.528895    0.568045    0.925597    0.071036    0.087129    0.020218    0.832620    0.778157    0.870012
2   0.978618    0.799159    0.461479    0.780529    0.118274    0.639921    0.143353    0.944669    0.521848    0.414662
3   0.264556    0.774234    0.456150    0.568434    0.018790    0.617635    0.612096    0.616934    0.943748    0.681820
4   0.359508    0.437032    0.697631    0.060225    0.666767    0.670638    0.210383    0.128926    0.315428    0.363711

Маскировка нескольких записей

A = A_orig.copy()
A[0, 0] = np.NAN
A[3, 1] = np.NAN
A[6, 3] = np.NAN

A_df = pd.DataFrame(A)
print A_df.head()


0   1   2   3   4   5   6   7   8   9
0   NaN 0.715189    0.602763    0.544883    0.423655    0.645894    0.437587    0.891773    0.963663    0.383442
1   0.791725    0.528895    0.568045    0.925597    0.071036    0.087129    0.020218    0.832620    0.778157    0.870012
2   0.978618    0.799159    0.461479    0.780529    0.118274    0.639921    0.143353    0.944669    0.521848    0.414662
3   0.264556    NaN 0.456150    0.568434    0.018790    0.617635    0.612096    0.616934    0.943748    0.681820
4   0.359508    0.437032    0.697631    0.060225    0.666767    0.670638    0.210383    0.128926    0.315428    0.363711

Определение матриц W и H

K = 4
W = np.abs(np.random.uniform(low=0, high=1, size=(M, K)))
H = np.abs(np.random.uniform(low=0, high=1, size=(K, N)))
W = np.divide(W, K*W.max())
H = np.divide(H, K*H.max())

pd.DataFrame(W).head()

0   1   2   3
0   0.078709    0.175784    0.095359    0.045339
1   0.006230    0.016976    0.171505    0.114531
2   0.135453    0.226355    0.250000    0.054753
3   0.167387    0.066473    0.005213    0.191444
4   0.080785    0.096801    0.148514    0.209789

pd.DataFrame(H).head()

0   1   2   3   4   5   6   7   8   9
0   0.074611    0.216164    0.157328    0.003370    0.088415    0.037721    0.250000    0.121806    0.126649    0.162827
1   0.093851    0.034858    0.209333    0.048340    0.130195    0.057117    0.024914    0.219537    0.247731    0.244654
2   0.230833    0.197093    0.084828    0.020651    0.103694    0.059133    0.033735    0.013604    0.184756    0.002910
3   0.196210    0.037417    0.020248    0.022815    0.171121    0.062477    0.107081    0.141921    0.219119    0.185125

Определение стоимости, которую мы хотим минимизировать

def cost(A, W, H):
    from numpy import linalg
    WH = np.dot(W, H)
    A_WH = A-WH
    return linalg.norm(A_WH, 'fro')

Однако, поскольку в A отсутствуют элементы, мы должны определить стоимость в терминах элементов, присутствующих в A.

def cost(A, W, H):
    from numpy import linalg
    mask = pd.DataFrame(A).notnull().values
    WH = np.dot(W, H)
    WH_mask = WH[mask]
    A_mask = A[mask]
    A_WH_mask = A_mask-WH_mask
    # Since now A_WH_mask is a vector, we use L2 instead of Frobenius norm for matrix
    return linalg.norm(A_WH_mask, 2)

Давайте просто попробуем увидеть стоимость начального набора значений W и H, которые мы случайно присвоили.

cost(A, W, H)

7.3719938519859509

Альтернативная процедура NNLS

num_iter = 1000
num_display_cost = max(int(num_iter/10), 1)
from scipy.optimize import nnls

for i in range(num_iter):
    if i%2 ==0:
        # Learn H, given A and W
        for j in range(N):
            mask_rows = pd.Series(A[:,j]).notnull()
            H[:,j] = nnls(W[mask_rows], A[:,j][mask_rows])[0]
    else:
        for j in range(M):
            mask_rows = pd.Series(A[j,:]).notnull()
            W[j,:] = nnls(H.transpose()[mask_rows], A[j,:][mask_rows])[0]
    WH = np.dot(W, H)
    c = cost(A, W, H)
    if i%num_display_cost==0:
        print i, c

0 4.03939072472
100 2.38059096458
200 2.35814781954
300 2.35717011529
400 2.35711130357
500 2.3571079918
600 2.35710729854
700 2.35710713129
800 2.35710709085
900 2.35710708109



A_pred = pd.DataFrame(np.dot(W, H))
A_pred.head()

0   1   2   3   4   5   6   7   8   9
0   0.564235    0.677712    0.558999    0.631337    0.536069    0.621925    0.629133    0.656010    0.839802    0.545072
1   0.788734    0.539729    0.517534    1.041272    0.119894    0.448402    0.172808    0.658696    0.493093    0.825311
2   0.749886    0.575154    0.558981    0.931156    0.270149    0.502035    0.287008    0.656178    0.588916    0.741519
3   0.377419    0.743081    0.370408    0.637094    0.071684    0.529433    0.767696    0.628507    0.832910    0.605742
4   0.458661    0.327143    0.610012    0.233134    0.685559    0.377750    0.281483    0.269960    0.468756    0.114950

Давайте просмотрим значения замаскированных записей.

A_pred.values[~pd.DataFrame(A).notnull().values]

array([ 0.56423481,  0.74308143,  0.10283106])

Original values were:

A_orig[~pd.DataFrame(A).notnull().values]

array([ 0.5488135 ,  0.77423369,  0.13818295])