矩阵在单纯形上的投影

Question

我在理解这段基于输出的代码时遇到了问题，我猜它会计算矩阵的特征向量

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）