There are bugs in 'auction_lap.py' but now I have fixed them

[563925743]

Hi Ben,

Thanks for sharing about the python implementation of the auction algorithm! I found some bugs in 'auction_lap.py' that actually make the codes output wrong results. I fixed them by myself and attached my codes. I am new to Pytorch. Maybe you can make my revised version better:)

1. In line 67, (curr_ass.view(-1, 1) == have_bidder.view(1, -1)).sum(dim=1) actually outputs an array containing 0 and 1 instead of the indices of assigned bids. At each round, it actually makes the first two elements of curr_ass be -1. It can be fixed by adding .nonzero() to the tail: (curr_ass.view(-1, 1) == have_bidder.view(1, -1)).sum(dim=1).nonzero()

2. The case when we only have 1 unassigned bidder at some round needs to be considered. Otherwise, we will get some dimensional errors. In my codes, I simply deal with this case separately.

I also tested my codes by solving a 2D points matching problem via auction algorithm and LP optimization using cvxpy package. They output the same results.

Here is the revised auction_lap.py:

#!/usr/bin/env python

from __future__ import division, print_function

import sys

import torch

def auction_lap(X, eps, compute_score=True):
    """
        X: n-by-n matrix w/ integer entries
        eps: "bid size" -- smaller values means higher accuracy w/ longer runtime
    """
   
    X = torch.from_numpy(X).float()
    # --
    # Init
    
    cost     = torch.zeros((1, X.shape[1]))
    curr_ass = torch.zeros(X.shape[0]).long() - 1
    bids     = torch.zeros(X.shape)
    
    # if X.is_cuda:
    #     cost, curr_ass, bids = cost.cuda(), curr_ass.cuda(), bids.cuda()
    
    
    counter = 0
    while (curr_ass == -1).any():
        counter += 1
        
        # --
        # Bidding
        
        unassigned = (curr_ass == -1).nonzero().squeeze()
        
        value = X[unassigned] - cost
        top_value, top_idx = value.topk(2, dim=1)
        
        first_idx = top_idx[:,0]
        first_value, second_value = top_value[:,0], top_value[:,1]
        
        bid_increments = first_value - second_value + eps
        
        bids_ = bids[unassigned]
        bids_.zero_()
        if unassigned.dim()==0:
            high_bids, high_bidders = bid_increments, unassigned
            cost[:,first_idx] += high_bids
            curr_ass[(curr_ass == first_idx).nonzero()] = -1
            curr_ass[high_bidders] = first_idx
        else:
            bids_.scatter_(
                dim=1,
                index=first_idx.contiguous().view(-1, 1),
                src=bid_increments.view(-1, 1)
            )
        
        
            have_bidder = (bids_ > 0).int().sum(dim=0).nonzero()
            
            high_bids, high_bidders = bids_[:,have_bidder].max(dim=0)
            
            high_bidders = unassigned[high_bidders.squeeze()]
            
            cost[:,have_bidder] += high_bids
            
            # curr_ass[(curr_ass.view(-1, 1) == have_bidder.view(1, -1)).sum(dim=1)] = -1
            ind = (curr_ass.view(-1, 1) == have_bidder.view(1, -1)).sum(dim=1).nonzero()
            curr_ass[ind] = -1
            
            curr_ass[high_bidders] = have_bidder.squeeze()
        
 
    score = None
    if compute_score:
        score = int(X.gather(dim=1, index=curr_ass.view(-1, 1)).sum())
    
    return score, curr_ass.numpy(), counter

Here is my testing code:

import timeit

import cvxpy as cp
import numpy as np
from scipy.sparse import lil_matrix


from auction_lap import auction_lap


np.random.seed(0)
n = 20
P = np.random.rand(n,2)
Q = np.random.rand(n,2)

start = timeit.default_timer()
C = lil_matrix((n,n))
for k in range(n):
    C[k,:] = -np.linalg.norm(P[k,:] - Q,axis=1)**2
    
C = C.toarray()
start = timeit.default_timer()
_, curr_ass, counter = auction_lap(C, eps=1e-6, compute_score=False)
stop = timeit.default_timer()
print('Autction time: ',stop-start)
M = np.array([np.linspace(0,n-1,n).astype(int).tolist(),curr_ass.tolist()])
S_auc = lil_matrix((n,n))
S_auc[M[1,:],M[0,:]] = 1
S_auc = S_auc.toarray()

start = timeit.default_timer()

S = cp.Variable((n,n))
constraints = [S >=0 , S @ np.ones(n) == 1, S.T @ np.ones(n) == 1]
prob = cp.Problem(
            cp.Maximize(cp.trace(Q.T @ S @ P)), 
            constraints
        )
prob.solve()

stop = timeit.default_timer() 
print('LP time: ',stop-start)
   
sol = S.value
sol[sol<1e-3] = 0
sol[sol>0.999] = 1

print('Difference between S_auc and sol: ',np.sum(np.abs(S_auc-sol)))

stop = timeit.default_timer()

Best,
Fengyu