Python linsys: Add normalization

[bamos]

Depends on: cvxgrp/scs#112

Example Usage

def normalize_a_cb(boundaries, scale):
    global A

    D_all = np.ones(ncon_cone)
    E_all = np.ones(nz_cone)

    min_scale, max_scale = (1e-4, 1e4)
    n_passes = 10

    for i in range(n_passes):
        D = np.sqrt(sla.norm(A, float('inf'), axis=1))
        E = np.sqrt(sla.norm(A, float('inf'), axis=0))
        D[D < min_scale] = 1.0
        E[E < min_scale] = 1.0
        D[D > max_scale] = max_scale
        E[E > max_scale] = max_scale
        start = boundaries[0]
        for delta in boundaries[1:]:
            D[start:start+delta] = D[start:start+delta].mean()
            start += delta
        A = sp.diags(1/D).dot(A).dot(sp.diags(1/E))
        D_all *= D
        E_all *= E

    mean_row_norm = sla.norm(A, 2, axis=1).mean()
    mean_col_norm = sla.norm(A, 2, axis=0).mean()
    A *= scale

    return D_all, E_all, mean_row_norm, mean_col_norm

def un_normalize_a_cb(D, E):
    global A
    A = sp.diags(D).dot(A).dot(sp.diags(E))

Output compared to vanilla SCS

D[0] = 1.000000, D[1] = 1.403820, D[2] = 0.594797, D[3] = 0.782981, D[4] = 0.
782981, D[5] = 0.782981, D[6] = 0.782981,
SCS(norm) D = 2.403483

E[0] = 1.209977, E[1] = 1.382466, E[2] = 2.373379,
SCS(norm) E = 3.001363
Setup time: 5.24e-04s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 2.50e+00  1.36e+00  8.00e-01 -2.00e+00  2.00e+00  8.73e-16  1.78e-04
     1| 4.64e-01  1.15e-01  2.73e-01  1.06e+00  3.92e-01  1.57e-16  6.93e-02
     2| 9.64e-02  2.53e-02  1.90e-03  5.09e-01  5.13e-01  9.87e-18  6.99e-02
     3| 5.11e-02  1.37e-02  5.33e-02  4.57e-01  5.65e-01  4.30e-17  7.00e-02
     4| 4.53e-02  1.48e-02  5.18e-02  4.65e-01  5.70e-01  3.52e-17  7.00e-02
     5| 3.38e-02  5.35e-03  6.62e-04  5.56e-01  5.54e-01  2.04e-16  7.00e-02
     6| 9.97e-03  5.62e-03  5.55e-03  5.56e-01  5.45e-01  1.14e-16  7.01e-02
     7| 5.78e-03  4.63e-03  4.43e-03  5.53e-01  5.44e-01  6.56e-17  7.01e-02
     8| 5.83e-03  2.28e-03  1.11e-03  5.51e-01  5.49e-01  8.34e-17  7.01e-02
     9| 1.89e-03  7.75e-04  1.00e-03  5.57e-01  5.55e-01  1.70e-16  7.02e-02
    10| 1.89e-03  7.75e-04  1.00e-03  5.57e-01  5.55e-01  1.70e-16  7.02e-02
----------------------------------------------------------------------------



D[0] = 1.000000, D[1] = 1.403820, D[2] = 0.594797, D[3] = 0.782981, D[4] = 0.
782981, D[5] = 0.782981, D[6] = 0.782981,
SCS(norm) D = 2.403483

E[0] = 1.209977, E[1] = 1.382466, E[2] = 2.373379,
SCS(norm) E = 3.001363
Setup time: 2.99e-02s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 2.50e+00  1.36e+00  8.00e-01 -2.00e+00  2.00e+00  8.73e-16  5.69e-04
     1| 4.64e-01  1.15e-01  2.73e-01  1.06e+00  3.92e-01  3.48e-17  8.37e-04
     2| 9.64e-02  2.53e-02  1.90e-03  5.09e-01  5.13e-01  7.20e-17  1.14e-03
     3| 5.11e-02  1.37e-02  5.33e-02  4.57e-01  5.65e-01  8.36e-17  1.40e-03
     4| 4.53e-02  1.48e-02  5.18e-02  4.65e-01  5.70e-01  8.23e-17  1.70e-03
     5| 3.38e-02  5.35e-03  6.62e-04  5.56e-01  5.54e-01  5.68e-17  1.95e-03
     6| 9.97e-03  5.62e-03  5.55e-03  5.56e-01  5.45e-01  6.00e-17  2.25e-03
     7| 5.78e-03  4.63e-03  4.43e-03  5.53e-01  5.44e-01  6.20e-17  2.59e-03
     8| 5.83e-03  2.28e-03  1.11e-03  5.51e-01  5.49e-01  1.57e-17  2.84e-03
     9| 1.89e-03  7.75e-04  1.00e-03  5.57e-01  5.55e-01  1.52e-17  3.08e-03
    10| 1.89e-03  7.75e-04  1.00e-03  5.57e-01  5.55e-01  1.52e-17  3.33e-03
----------------------------------------------------------------------------

Memory Usage

I ran the same tests from #11.
The version with Python callbacks is using more memory since there are Python objects moving around. It doesn't seem like there are significant memory leaks because there are some iterations where the memory consumption doesn't increase at all.

Setup

data, chain, inv_data = prob.get_problem_data('SCS')
print('\n\n--- vanilla scs, same data ---')
gc.collect(); print('  + {} bytes'.format(process.memory_info().rss))
for _ in range(10):
    sol = scs.solve(
        data, cones, verbose=False,
        use_indirect=False, normalize=True,
        max_iters=10,
        scale=0.5,
    )
    gc.collect(); print('  + {} bytes'.format(process.memory_info().rss))


print('\n\n--- vanilla scs, random data ---')
init_m = process.memory_info().rss
gc.collect(); print('  + {} bytes'.format(init_m))
for _ in range(20):
    _G.value = npr.randn(nineq, nz)
    data, chain, inv_data = prob.get_problem_data('SCS')
    sol = scs.solve(
        data, cones, verbose=False,
        use_indirect=False, normalize=True,
        max_iters=10,
        scale=0.5,
    )
    gc.collect()
    m = process.memory_info().rss
    print('  + {} bytes (+ {} bytes)'.format(m, m-init_m))
    init_m = m

Output

--- vanilla scs, same data ---
  + 108650496 bytes
  + 112267264 bytes
  + 112267264 bytes
  + 112267264 bytes
  + 112267264 bytes
  + 112267264 bytes
  + 112267264 bytes
  + 112267264 bytes
  + 112267264 bytes
  + 112267264 bytes
  + 112267264 bytes


--- vanilla scs, random data ---
  + 112267264 bytes
  + 112295936 bytes (+ 28672 bytes)
  + 112320512 bytes (+ 24576 bytes)
  + 112349184 bytes (+ 28672 bytes)
  + 112398336 bytes (+ 49152 bytes)
  + 112431104 bytes (+ 32768 bytes)
  + 112455680 bytes (+ 24576 bytes)
  + 112504832 bytes (+ 49152 bytes)
  + 112537600 bytes (+ 32768 bytes)
  + 112541696 bytes (+ 4096 bytes)
  + 112545792 bytes (+ 4096 bytes)
  + 112545792 bytes (+ 0 bytes)
  + 112545792 bytes (+ 0 bytes)
  + 112545792 bytes (+ 0 bytes)
  + 112562176 bytes (+ 16384 bytes)
  + 112562176 bytes (+ 0 bytes)
  + 112562176 bytes (+ 0 bytes)
  + 112566272 bytes (+ 4096 bytes)
  + 112566272 bytes (+ 0 bytes)
  + 112566272 bytes (+ 0 bytes)
  + 112578560 bytes (+ 12288 bytes)


--- scs with python cbs, same data ---                               [3/1942]
  + 112652288 bytes
/Users/bamos/anaconda3/lib/python3.6/site-packages/scipy/sparse/linalg/dsolv$
/linsolve.py:296: SparseEfficiencyWarning: splu requires CSC matrix format
  warn('splu requires CSC matrix format', SparseEfficiencyWarning)
  + 113025024 bytes
  + 113029120 bytes
  + 113033216 bytes
  + 113037312 bytes
  + 113041408 bytes
  + 113045504 bytes
  + 113049600 bytes
  + 113049600 bytes
  + 113053696 bytes
  + 113057792 bytes


--- scs with python cbs, random data ---
  + 113057792 bytes
  + 113115136 bytes (+ 57344 bytes)
  + 113180672 bytes (+ 65536 bytes)
Error in AA type 1, iter: 9, info: 0, norm 1.40e+04
  + 113233920 bytes (+ 53248 bytes)
  + 113274880 bytes (+ 40960 bytes)
  + 113315840 bytes (+ 40960 bytes)
  + 113340416 bytes (+ 24576 bytes)
  + 113352704 bytes (+ 12288 bytes)
  + 113360896 bytes (+ 8192 bytes)
  + 113393664 bytes (+ 32768 bytes)
  + 113397760 bytes (+ 4096 bytes)
  + 113426432 bytes (+ 28672 bytes)
  + 113430528 bytes (+ 4096 bytes)
  + 113438720 bytes (+ 8192 bytes)
  + 113455104 bytes (+ 16384 bytes)
  + 113459200 bytes (+ 4096 bytes)
  + 113463296 bytes (+ 4096 bytes)
  + 113463296 bytes (+ 0 bytes)
  + 113467392 bytes (+ 4096 bytes)
  + 113471488 bytes (+ 4096 bytes)
  + 113475584 bytes (+ 4096 bytes)

[bodono]

LGTM. You can update the scs submodule to pull in your latest change there as part of this PR, once that is done let me know and I'll merge this in.

[bamos]

Added