Scalings of Matricies Satisfying Line-Product Constraints (Python)

line_product_scaling.py

# Python implementation of algorithms described in
# "Scalings of Matrices Satisfying Line-Product Constraints and Generalizations" (1992)
# Paper by Uriel G. Rothblum, code by @tjbanks and Tung Nguyen
# https://core.ac.uk/download/pdf/82705928.pdf
# Implement 2D matrix completion with input as U,V, A_prime and output as A_exp

import numpy as np

def completion(U = None,V = None,A_prime = None): 
    # Situate unknown entries as 1.
    A_prime[A_prime == 0] = 1.0
    # Scale back to the exponential space
    A_exp = np.multiply(A_prime,1.0) / (U * np.transpose(V))
    return A_exp

# Implement A_prime to potentially put into another architecture. Input is A, and outputs are
# U, V, and A_prime

def get_A_prime(A = None): 
    ##### Step 0 in the paper
    m,n = A.shape
    # Initiate sigma(A)
    rho_sign = np.ones(A.shape)
    # Initiate sigma_row and sigma_col to get the number of nonzeros inside each row and column
    sigma_row = np.zeros([m,1])
    sigma_col = np.zeros([n,1])
    # Setup sigma(A) that indicates known and unknown entries.
    rho_sign = A != 0

    # Get the number of nonzeros inside each row and column
    for i in np.arange(0,m).reshape(-1):
        sigma_row[i] = sum(rho_sign[i,:])

    for j in np.arange(0,n).reshape(-1):
        sigma_col[j] = sum(rho_sign[:,j])

    # Take logarithm of matrix
    a = np.log(A)
    
    # After log, all 0 values will be -inf, so we set them to 0
    a[a == -np.inf] = 0.0

    # Initiate U and V
    u = s = np.zeros([m,1])
    v = t = np.zeros([n,1])
    # Initiate the stopping point for both step 1 and 2 in the paper.
    epsilon = 10 ** - 22
    sumv_row = sumv_col = 0.0
    prev_row = prev_col = 1.0
    ##### Step 1 and 2 in the paper
    trial = 0
    while np.abs(prev_row - sumv_row) > epsilon or np.abs(prev_col - sumv_col) > epsilon:

        # Starting the first iterative step
        for col in np.arange(0,n).reshape(-1):
            # Get the sum by columns first
            sig_size = sigma_col[col][0]
            #update rho_col
            if sig_size > 0.0:
                rho_col = - sum(a[:,col]) / sig_size
                # update a_ij
                a[:,col] += rho_col * rho_sign[:,col]
                v[col] += rho_col
        for row in np.arange(0,m).reshape(-1):
            # Get the sum by rows first
            sig_size = sigma_row[row][0]
            if sig_size > 0.0:
                #update rho_row
                rho_row = - sum(a[row,:]) / sig_size
                a[row,:] += rho_row * rho_sign[row,:]
                u[row] += rho_row
        prev_row = sumv_row
        prev_col = sumv_col
        # Sum is added to get only (1 x column) remain
        sumv_row = np.var(sum(a))
        sumv_col = np.var(sum(np.transpose(a)))
        trial += 1
        # print('Iteration %d\n' % (trial))


    ##### Method 1: Go to exp space and multiply by inverse
    # Finally, we have all the predicted elements inside matrix a.

    # Avoid blowing up coefficients
    A_prime = np.multiply(np.exp(a),np.sign(A))
    U = np.exp(u)
    V = np.exp(v)
    return U,V,A_prime

a = np.array([[.5,.3,.8],[.1,.9,.8]])
U,V,A_prime = get_A_prime(a)
print('A: {}\n'.format(a))
print('U: {}\n'.format(U))
print('V: {}\n'.format(V))
print('A_prime: {}\n'.format(A_prime))

Example run:

>python line_product_scaling.py

A: [[0.5 0.3 0.8]
 [0.1 0.9 0.8]]

U: [[0.9183859 ]
 [1.08886689]]

V: [[4.47213595]
 [1.9245009 ]
 [1.25      ]]

A_prime: [[2.05357331 0.53023035 0.9183859 ]
 [0.48695608 1.88597277 1.08886689]]