Проекция матрицы на симплекс - PullRequest
Новая
       49

Проекция матрицы на симплекс

0 голосов
/ 20 января 2020

У меня проблема с пониманием этого фрагмента кода, который основан на выводе, я думаю, он вычисляет собственный вектор матрицы. 

def simplexProj(y):
    """
    Given y,  computes its projection x* onto the simplex

          Delta = { x | x >= 0 and sum(x) <= 1 },

    that is, x* = argmin_x ||x-y||_2  such that  x in Delta.

    x = SimplexProj(y)

    ****** Input ******
    y    : input vector.

    ****** Output ******
    x    : projection of y onto Delta.
    """

    if len(y.shape) == 1:  # Reshape to (1,-1) if y is a vector.
        y = y.reshape(1, -1) # row vector

    x = y.copy()
    x[x < 0] = 0 #element within the matrix that is negative will be replaced with 0, python2 feature
    K = np.flatnonzero(np.sum(x, 0) > 1) #return indices that are non-zero in the flattened version of a ; sum of each column
    # K gives the column index for column that has colum sum>1, True = 1, False = 0
    x[:, K] = blockSimplexProj(y[:, K])
    return x

def blockSimplexProj(y):
    """ Same as function SimplexProj except that sum(max(Y,0)) > 1. """

    r, c = y.shape
    ys = -np.sort(-y, axis=0) #sort each column of the matrix with biggest entry on the first row
    mu = np.zeros(c, dtype=float)
    S = np.zeros((r, c), dtype=float)

    for i in range(1, r): #1st to r-1th row
        S[i, :] = np.sum(ys[:i, :] - ys[i, :], 0)
        print(S)
        colInd_ge1 = np.flatnonzero(S[i, :] >= 1)
        colInd_lt1 = np.flatnonzero(S[i, :] < 1)
        if len(colInd_ge1) > 0:
            mu[colInd_ge1] = (1 - S[i - 1, colInd_ge1]) / i - ys[i - 1, colInd_ge1]
        if i == r:
            mu[colInd_lt1] = (1 - S[r, colInd_lt1]) / (r + 1) - ys[r, colInd_lt1]

    x = y + mu
    x[x < 0] = 0
    return x

Я немного озадачен пошаговым вычислением матрицы S потому что согласно коду строка первой строки S должна быть равна 0. Возьмем, к примеру, матрицу A = np.array([[25,70,39,10,80],[12,45,32,89,43],[67,24,84,39,21],[0.1,0.2,0.3,0.035,0.06]]) 3 итерации (i=1,2,3) вычисляются, как и ожидалось, но затем есть дополнительный шаг, который похоже, возвращает S в качестве основы собственных векторов. Было бы здорово, если бы кто-нибудь помог мне разобраться в этой проблеме. Также я # не уверен, как называется этот алгоритм (как вычисляется S)

...