Optimized Resolvent Splitting Examples

This notebook generates the figures in our paper, which outlines the use of a semi-definite programming framework for optimizing splitting algorithms which solve the following monotone inclusion problem:

\begin{equation} 0 \in \sum_{i=1}^{n} A_{i}(x). \end{equation}

It does so by finding valid matrices $L$ and $W$ for the following set of iterations, where bolded entries are lifted, i.e. $\mathrm{x} = (x_1, \dots, x_n)$:

\begin{align} \mathrm{x} &= J_{\mathrm{A}}\left(\mathrm{v}^{k} + \mathrm{L} \mathrm{x}\right)\\ \mathrm{v}^{k+1} &= \mathrm{v}^{k} - \gamma \mathrm{W} \mathrm{x}. \end{align}

$L$ is lower diagonal, so we work with a symmetrized matrix $Z = 2I - L - L^T$

import numpy as np
import oars
from oars.matrices import *
from oars.utils import *
from oars.pep import *
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
# matplotlib.use('nbagg')
import seaborn as sns
sns.set()
from matplotlib.colors import CenteredNorm
import cvxpy as cvx
from collections import defaultdict

Our first example builds the matrices, which can be thought of as the graph Laplacians of communication graphs for the nodes in a distributed setting, to avoid specific communication paths.

We suppose a set of two clusters (0-2) and (3-5) which can only communicate within clusters and via a link between 0-3.

We then plot a Gantt chart for the distributed execution of the resulting algorithm.

We take as our objective function here the spectral norm of $Z-W$, which promotes graph similarity.

# Data for matrices in cluster example (6x6)
n = 6
Z_fixed = {(r,c):0 for r in range(4,6) for c in range(3)}
Z_fixed[(3,1)] = 0
Z_fixed[(3,2)] = 0

Z, W = getMinSpectralDifference(n, fixed_Z=Z_fixed, fixed_W=Z_fixed)

# Figure 3
t = 0.016*np.ones(n)
l = 0.00025*np.ones((n, n))
l[0, 3] = 0.01
l[3, 0] = 0.01  
fig = getGantt(t, l, Z, W, title="Cluster Algorithm Parallelism", itrs=5)
fig.show()
#if you want to save the figure as a separate file
saving = False # set to True to save the figures
if saving:
    import os
    if ~os.path.exists('figs'): os.makedirs('figs')
    fig.write_image("figs/fig_gantt.pdf")

In this example, we visually depict the different matrices generates by changing the objective function between maximizing algebraic connectivity, minimizing SLEM, minimizing total effective resistance, and minimizing the spectral norm of the difference of the matrices.

We also demonstrate the use of the getBlockMin function to build the constraint set for the SDP to guarantee that the nodes can execute in blocks. If the compute and communication times are constant and uniform, this guarantees that the algorithm can execute with the highest possible degree of parallelism possible for this type of splitting.

# Figure 4 helper
def cplot(W,Z):
    vmin = np.floor(np.min(np.minimum(W, Z)))
    cm = plt.cm.coolwarm
    fig, ax = plt.subplots(1, 2, figsize=(12, 6))
    x = -np.arange(n)
    y = np.arange(n)
    ax[0].pcolormesh(x, y, W, norm=CenteredNorm(), cmap=cm) #, vmin=vmin, vmax=2)
    ax[0].set_title(r"$W$")
    ax[0].axis('off')
    ax[0].set_aspect('equal')

    pc = ax[1].pcolormesh(x, y, Z, norm=CenteredNorm(), cmap=cm)
    ax[1].set_title(r"$Z$")
    ax[1].axis('off')
    ax[1].set_aspect('equal')

    # Add colorbar
    plt.tight_layout()
    
    fig.colorbar(pc, ax=ax.ravel().tolist())
    return fig

# Figure 4
n = 10
m = 2
Z, W = getBlockMin(n, m, builder=getMaxConnectivity)

fig = cplot(W, Z)

# Save as pdf
if saving: fig.savefig("figs/fig_block_fiedler.pdf")

# Figure 4
Z, W = getBlockMin(n, m, builder=getMinSLEM)

fig = cplot(W, Z)
if saving: fig.savefig("figs/fig_block_slem.pdf")

# Figure 4
Z, W = getBlockMin(n, m, builder=getMinResist)

fig = cplot(W, Z)
if saving: fig.savefig("figs/fig_block_resist.pdf")

# Figure 4
Z, W = getBlockMin(n, m, builder=getMinSpectralDifference)

fig = cplot(W, Z)
if saving: fig.savefig("figs/fig_block_lw.pdf")

This example shows a comparison of the execution of the block designs when the computation times are not uniform, and the use of the getMinIteration function to generate designs which minimize the iteration time.

# Figure 5 (hangs if executed in Windows using Kaleido 0.1.0post1)
# try pip install kaleido==0.1.0post1
n = 6
t = np.ones(n)
l = np.ones((n,n))/10
t[0] = 6
t[2] = 6
t[3] = 6

for i in range(1, n//2+1):
    Zb, Wb = getBlockMin(n, i, builder=getMinResist)
    if i == 1:
        itrs = 6
    else:
        itrs = 7
    fig = getGantt(t, l, Zb, Wb, title=str(n//i)+"-Block", itrs=itrs)
    
    if saving: fig.write_image("figs/fig_misdp"+str(n//i)+"block.pdf")
    fig.show()

Z, W = getMinIteration(n, t=t, l=l, builder=getMinResist)

fig = getGantt(t, l, Z, W, title="Minimum Iteration", itrs=12)

if saving: fig.write_image("figs/fig_misdp_min_cycle.pdf")
fig.show()

In several of the following examples, we use the PEP formulation to find the contraction factor for the algorithm designs in two scenarios:

1. $(n)$ 2-Lipschitz 1-strongly monotone operators

2. $(n-1)$ 2-Lipschitz 1-strongly monotone operators, and 1 unrestricted maximal monotone operator (this simulates having the subdifferential of an indicator function over some convex set)

In the figure 6 example, we use a fixed step size of 0.5.

We then compare the different algorithm designs, including the spectral objectives using our SDP, and some existing designs in the literature.

from warnings import warn
if 'MOSEK' not in cvx.installed_solvers():
    warn('''These experiments were conducted with Mosek as the underlying conic solver for CVXPY.  Mosek was not detected on your system. Running without Mosek will generate different results than those in the paper, including lack of convergence on larger problems.''')

# Figure 6 - contraction rate results
# fig_smoothstronglyconvex.ipynb
# Function to generate list of all test matrices
def getTestMatrices(n):
    cases = []
    Ws = []
    names = ['Full', 'MT']
    Lf, Wf = getFull(n)
    Mf = getMfromWCholesky(Wf)
    cases.append((Lf, Mf))
    Ws.append(Wf)
    Lmt, Wmt = getMT(n)
    Mmt = getMfromWCholesky(Wmt)
    cases.append((Lmt, Mmt))
    Ws.append(Wmt)
    i = n//2
    for obj in [getMinSLEM, getMinResist, getMaxConnectivity]:
        L, W = getBlockMin(n, i, builder=obj)
        M = getMfromWCholesky(W)
        cases.append((L, M))
        Ws.append(W)
        title = f"2-Block "+ obj.__name__
        names.append(title)
        
    return cases, names, Ws
    
def getMatrixContraction(cases, names, Lc=2, mu=1, gamma=0.5):
    '''Get contraction factors for matrix cases

    Args:
        cases (list): list of cases with L and M
        names (list): list of names
        Lc (float): Lipschitz constant
        mu (float): strong convexity parameter
        gamma (float): step size

    Returns:
        taus (dict): dictionary of contraction factors
    '''

    taus = {}
    for i, (L, M) in enumerate(cases):
        n = L.shape[0]
        
        ls = np.ones(n)*Lc
        mus = np.ones(n)*mu
        
        tau = getReducedContractionFactor(L, M, ls, mus)
        taus[names[i]] = tau
        print(names[i], tau)
    return taus

# Figure 6 helper
def testFixedRate(start=4, end=31, lc=2, mc=1):
    tau_results = defaultdict(list) 
    for n in range(start, end, 2):
        print(f"Testing optimal rate on {n} nodes")
        ls = np.ones(n) * lc
        mus = np.ones(n) * mc
        cases, names, Ws = getTestMatrices(n)
        for i, (L, M) in enumerate(cases):
            tau = getReducedContractionFactor(L, M, ls, mus)
            tau_results[names[i]].append(tau)


    # Convert results to df
    tau_df = pd.DataFrame(tau_results)
    tau_df.index = range(start, end, 2)

    return tau_df

# Figure 6a
tau_df = testFixedRate(start=6, end=31)

Testing optimal rate on 6 nodes
Testing optimal rate on 8 nodes
Testing optimal rate on 10 nodes
Testing optimal rate on 12 nodes
Testing optimal rate on 14 nodes
Testing optimal rate on 16 nodes
Testing optimal rate on 18 nodes
Testing optimal rate on 20 nodes
Testing optimal rate on 22 nodes
Testing optimal rate on 24 nodes
Testing optimal rate on 26 nodes
Testing optimal rate on 28 nodes
Testing optimal rate on 30 nodes

# Plot figure 6a
fig, ax = plt.subplots()
tau_df.plot(ax=ax)

ax.set_xlabel("n")
ax.set_ylabel("Contraction")
ax.set_title(r"Fixed step size contraction with $n$ Lipschitz strongly monotone operators")
ax.legend(['Full', 'MT', '2-Block Min SLEM', r'2-Block Min Resistance (also min $\|Z-W\|$)', '2-Block Max Fiedler'], loc=(0.25,0.4))
if saving: fig.savefig("figs/fig6a.pdf", format="pdf", bbox_inches="tight")

# Figure 6b helper
def testFixedRateProj(start=4, end=31, lc=2, mc=1):
    tau_results = defaultdict(list) 
    for n in range(start, end, 2):
        print(f"Testing optimal rate on {n} nodes")
        ls = np.ones(n) * lc
        ls[n-1] = np.inf
        mus = np.ones(n) * mc
        mus[n-1] = 0
        cases, names, Ws = getTestMatrices(n)
        for i, (Z, M) in enumerate(cases):
            tau = getReducedContractionFactor(Z, M, ls, mus)
            tau_results[names[i]].append(tau)


    # Convert results to df
    tau_df = pd.DataFrame(tau_results)
    tau_df.index = range(start, end, 2)

    return tau_df

# Figure 6b
tau_df = testFixedRateProj(start=6, end=31)

Testing optimal rate on 6 nodes
Testing optimal rate on 8 nodes
Testing optimal rate on 10 nodes
Testing optimal rate on 12 nodes
Testing optimal rate on 14 nodes
Testing optimal rate on 16 nodes
Testing optimal rate on 18 nodes
Testing optimal rate on 20 nodes
Testing optimal rate on 22 nodes
Testing optimal rate on 24 nodes
Testing optimal rate on 26 nodes
Testing optimal rate on 28 nodes
Testing optimal rate on 30 nodes

# Plot figure 6b
fig, ax = plt.subplots()
tau_df.plot(ax=ax)

ax.set_xlabel("n")
ax.set_ylabel("Contraction")
ax.set_title("Fixed step size contraction with an unrestricted monotone operator")
ax.legend(['Full', 'MT', '2-Block Min SLEM', r'2-Block Min Resistance (also min $\|Z-W\|$)', '2-Block Max Fiedler'], loc="lower right")
if saving: fig.savefig("figs/fig6b.pdf", format="pdf", bbox_inches="tight")

In figure 7 we repeat the analysis, but using an optimal step size we derive from an extension of the dual of the PEP. This also demonstrates the use of the getReducedContractionOptGamma function, which generates a step size which minimizes the contraction factor for the reduced version of the algorithm:

\begin{align} \mathrm{x} &= J_{\mathrm{A}}\left(-\mathrm{M}^{T} \mathrm{z}^{k} + \mathrm{L} \mathrm{x}\right)\\ \mathrm{z}^{k+1} &= \mathrm{z}^{k} + \gamma \mathrm{M} \mathrm{x}. \end{align}

where $M^T M = W$

# Figure 7 helper
def testOptimalRate(start=4, end=31, lc=2, mc=1):
    tau_results = defaultdict(list) 
    for n in range(start, end, 2):
        print(f"Testing optimal rate on {n} nodes")
        ls = np.ones(n) * lc
        ms = np.ones(n) * mc
        cases, names, Ws = getTestMatrices(n)
        for i, (Z, M) in enumerate(cases):
            tau, gamma = getReducedContractionOptGamma(Z, M, ls, ms)
            tau_results[names[i]].append(tau)


    # Convert results to df
    tau_df = pd.DataFrame(tau_results)
    tau_df.index = range(start, end, 2)

    return tau_df

# Figure 7 helper
def plotTau(tau_df, title):
    fig, ax = plt.subplots()
    tau_df.plot(ax=ax)
    ax.set_xlabel("n")
    ax.set_ylabel("Contraction") #$\\tau$
    ax.set_title(title)
    ax.legend(['Full', 'MT', '2-Block Min SLEM', r'2-Block Min Resistance (also min $\|Z-W\|$)', '2-Block Max Fiedler'], loc="center right")
    #ax.legend(loc='center right')
    #plt.show()
    return fig

# Figure 7
tau_df = testOptimalRate(start=4, end=31)

if saving: tau_df.to_csv("figs/fig_opt_gamma_strong.csv", index=False)

Testing optimal rate on 4 nodes
Testing optimal rate on 6 nodes
Testing optimal rate on 8 nodes
Testing optimal rate on 10 nodes
Testing optimal rate on 12 nodes
Testing optimal rate on 14 nodes
Testing optimal rate on 16 nodes
Testing optimal rate on 18 nodes
Testing optimal rate on 20 nodes
Testing optimal rate on 22 nodes
Testing optimal rate on 24 nodes
Testing optimal rate on 26 nodes
Testing optimal rate on 28 nodes
Testing optimal rate on 30 nodes

# Figure 7
fig = plotTau(tau_df, "Contraction over $n$ Lipschitz strongly monotone operators")
if saving: fig.savefig("figs/fig_opt_gamma_strong.pdf", format="pdf", bbox_inches="tight")

# Figure 7 helper
def testOptimalRateProj(start=4, end=31, lc=2, mc=1):
    tau_results = defaultdict(list) 
    for n in range(start, end, 2):
        print(f"Testing optimal rate on {n} nodes")
        ls = np.ones(n) * lc
        ls[n-1] = np.inf
        mus = np.ones(n) * mc
        mus[n-1] = 0
        cases, names, Ws = getTestMatrices(n)
        for i, (Z, M) in enumerate(cases):
            tau, _ = getReducedContractionOptGamma(Z, M, ls, mus)
            tau_results[names[i]].append(tau)


    # Convert results to df
    tau_df = pd.DataFrame(tau_results)
    tau_df.index = range(start, end, 2)

    return tau_df

# Figure 7b
tau_df = testOptimalRateProj(start=4, end=31)

if saving: tau_df.to_csv("figs/fig_opt_gamma_proj.csv", index=False)

Testing optimal rate on 4 nodes
Testing optimal rate on 6 nodes
Testing optimal rate on 8 nodes
Testing optimal rate on 10 nodes
Testing optimal rate on 12 nodes
Testing optimal rate on 14 nodes
Testing optimal rate on 16 nodes
Testing optimal rate on 18 nodes
Testing optimal rate on 20 nodes
Testing optimal rate on 22 nodes
Testing optimal rate on 24 nodes
Testing optimal rate on 26 nodes
Testing optimal rate on 28 nodes
Testing optimal rate on 30 nodes

# Figure 7b
fig = plotTau(tau_df, "Contraction with an unrestricted monotone operator")

if saving: fig.savefig("figs/fig7b.pdf", format="pdf", bbox_inches="tight")

In figure 8 we examine the impact of operator ordering on the contraction factor. In this example (and previous ones), we illustrate the use of some of the prebuilt matrix functions.

getMT(n) returns the Malitsky-Tam algorithm design over $n$ operators

getRyu(n) returns the Tam extension of the Ryu algorithm design over $n$ operators

getFull(n) returns the fully connected algorithm with equal weights over $n$ operators

We also use a helper function getMfromWCholesky to return a valid $M \in \mathbf{R}^{n-1 \times n}$ such that $M^T M = W$

# Figure 8 - operator ordering
Lc=2
mu=1


taus = defaultdict(list)
for n in range(4, 31, 2):
    print(f"Testing optimal rate on {n} nodes")
    Lb, Wb = getBlockMin(n, n//2)
    Mb = getMfromWCholesky(Wb)
    Lmt, Wmt = getMT(n)
    Mmt = getMfromWCholesky(Wmt)
    Lr, Wr = getRyu(n)
    Mr = getMfromWCholesky(Wr)
    Lmax, Wmax = getFull(n)
    Mmax = getMfromWCholesky(Wmax)

    cases = [(Lb, Mb), (Lmt, Mmt), (Lr, Mr), (Lmax, Mmax)]
    names = ["2-Block", "MT", "Ryu", "Full"]
    ls = np.ones(n)*Lc
    ls[n-1] = np.inf
    mus = np.ones(n)*mu
    mus[n-1] = 0
    for i, (Z, M) in enumerate(cases):
        tau, _ = getReducedContractionOptGamma(Z, M, ls, mus)
        taus[names[i]].append(tau)

    ls = np.ones(n)*Lc
    ls[0] = np.inf
    mus = np.ones(n)*mu
    mus[0] = 0
    for i, (Z, M) in enumerate(cases):
        tau, _ = getReducedContractionOptGamma(Z, M, ls, mus)        
        taus[names[i]+" reverse"].append(tau)

Testing optimal rate on 4 nodes
Testing optimal rate on 6 nodes
Testing optimal rate on 8 nodes
Testing optimal rate on 10 nodes
Testing optimal rate on 12 nodes
Testing optimal rate on 14 nodes
Testing optimal rate on 16 nodes
Testing optimal rate on 18 nodes
Testing optimal rate on 20 nodes
Testing optimal rate on 22 nodes
Testing optimal rate on 24 nodes
Testing optimal rate on 26 nodes
Testing optimal rate on 28 nodes
Testing optimal rate on 30 nodes

# Figure 8 - operator ordering
df = pd.DataFrame(taus)
df.index = range(4, 31, 2)

# Plot the results
colors = ['blue', 'green', 'red', 'purple']
styles = ['b-', 'g-', 'k-', 'm-', 'y--', 'c--', 'r--', 'p--']
styles = ['b-', 'g-', 'k-', 'm-', 'b*-', 'g*-', 'k*-', 'm*-']

if saving: df.to_csv("figs/fig_function_ordering_opt.csv", index=False)
df.plot(style=styles)
plt.xlabel("Number of nodes")
plt.ylabel("Contraction")
plt.title("Operator Ordering and Splitting Performance")

if saving: plt.savefig("figs/fig_function_order_opt.pdf", format="pdf", bbox_inches="tight")
plt.show()

Figure 9 examines the impact of block count in the $d$-Block designs on the reduced contraction factor.

# Figure 9 - block size

# Function to generate list of all test matrices
def getTestMatricesBlocks(n):
    cases = []
    names = ['Full', 'MT']
    objs = ['Full', 'MT']
    blocksizes = ['', '']
    blocklist = ['', '']
    Lf, Wf = getFull(n)
    Mf = getMfromWCholesky(Wf)
    cases.append((Lf, Mf))
    Lmt, Wmt = getMT(n)
    Mmt = getMfromWCholesky(Wmt)
    cases.append((Lmt, Mmt))
    objstring = ['MinSLEM', 'MinResist', 'getMaxConnectivity', 'Min |Z-W|']
    for i in range(2, n//2+1):
        if n % i != 0:
            continue
        for j, obj in enumerate([getMinSLEM, getMinResist, getMaxConnectivity, getMinSpectralDifference]):
            L, W = getBlockMin(n, i, builder=obj)
            M = getMfromWCholesky(W)
            cases.append((L, M))
            objs.append(objstring[j])
            blocksizes.append(i)
            blocklist.append(n//i)
            title = f"{n//i}-Block "+objstring[j]
            names.append(title)
        
    return cases, names, objs, blocksizes, blocklist

# Test function
def testMatrices(cases, names, lc=2, mu=1):
        
    taus = {}
    for i, (L, M) in enumerate(cases):
        ls = np.ones(len(L))*lc
        mus = np.ones(len(L))*mu
        
        tau, _ = getReducedContractionOptGamma(L, M, ls, mus)
        print(names[i], tau)
        taus[names[i]] = tau

    return taus

    

# Figure 9
results = []
for n in range(4, 25): #t < 120:
    print(n)
    cases, names, objs, blocksizes, blocks = getTestMatricesBlocks(n)
    taus = testMatrices(cases, names)
    # Each key in taus is a column, the row is the value of n
    for i in range(len(names)):
        k = names[i]
        v = taus[k]
        results.append({"n": n, "name": k, "obj": objs[i], "blocksize": blocksizes[i], "blocks":blocks[i], "tau": v})

df = pd.DataFrame(results)

if saving: df.to_csv("figs/taus_plus.csv", index=False)
df.set_index('n', inplace=True)

4
Full 0.42757264243870163
MT 0.8002281453197669
2-Block MinSLEM 0.5571563320360297
2-Block MinResist 0.42857143084394755
2-Block getMaxConnectivity 0.4723279099751601
2-Block Min |Z-W| 0.42857142974428214
5
Full 0.4476922209675467
MT 0.8613717840444194
6
Full 0.4592582042344319
MT 0.8973983347682897
3-Block MinSLEM 0.7135848361436024
3-Block MinResist 0.7182631442291106
3-Block getMaxConnectivity 0.7125769075299765
3-Block Min |Z-W| 0.7146237387649743
2-Block MinSLEM 0.5571562848352016
2-Block MinResist 0.4285714345044067
2-Block getMaxConnectivity 0.4884713576568497
2-Block Min |Z-W| 0.4285714294945957
7
Full 0.4667573414540386
MT 0.9210285463263077
8
Full 0.47200715470169285
MT 0.9373167138798192
4-Block MinSLEM 0.830902661865993
4-Block MinResist 0.8108953222170883
4-Block getMaxConnectivity 0.8002350026587245
4-Block Min |Z-W| 0.9414775041852568
2-Block MinSLEM 0.5571562864004252
2-Block MinResist 0.42857143242269075
2-Block getMaxConnectivity 0.4971297780170147
2-Block Min |Z-W| 0.42857157367446364
9
Full 0.47588407918368236
MT 0.9490542144444052
3-Block MinSLEM 0.7139323261605606
3-Block MinResist 0.7182631341824356
3-Block getMaxConnectivity 0.7138583593909309
3-Block Min |Z-W| 0.7154474793706733
10
Full 0.47886238699142825
MT 0.9578086304484189
5-Block MinSLEM 0.8909710612103936
5-Block MinResist 0.8675177048793344
5-Block getMaxConnectivity 0.869466506289027
5-Block Min |Z-W| 0.8841703860126198
2-Block MinSLEM 0.5571562868119783
2-Block MinResist 0.4285714375495388
2-Block getMaxConnectivity 0.5007530033325157
2-Block Min |Z-W| 0.4285714341245416
11
Full 0.48122091775576215
MT 0.9645090300428389
12
Full 0.4831342670605387
MT 0.9697458591678797
6-Block MinSLEM 0.9274863709906866
6-Block MinResist 0.9043223463177122
6-Block getMaxConnectivity 0.9144355840496884
6-Block Min |Z-W| 0.9723015042825072
4-Block MinSLEM 0.8306043397184802
4-Block MinResist 0.8108995469242436
4-Block getMaxConnectivity 0.8002302736642987
4-Block Min |Z-W| 0.9732200751628457
3-Block MinSLEM 0.7138823830882524
3-Block MinResist 0.7182631630830557
3-Block getMaxConnectivity 0.7152305558438894
3-Block Min |Z-W| 0.7154244052790969
2-Block MinSLEM 0.557156288107003
2-Block MinResist 0.4285714382735192
2-Block getMaxConnectivity 0.5048842320957917
2-Block Min |Z-W| 0.42857143414684185
13
Full 0.48471718030959554
MT 0.973913351855411
14
Full 0.4860481904482468
MT 0.9772819066072854
7-Block MinSLEM 0.9487402257486269
7-Block MinResist 0.9275962948707096
7-Block getMaxConnectivity 0.9415108622982323
7-Block Min |Z-W| 0.9379760461799428
2-Block MinSLEM 0.5571563188192495
2-Block MinResist 0.42857144512733913
2-Block getMaxConnectivity 0.5053543890781282
2-Block Min |Z-W| 0.4285714453701384
15
Full 0.4871828194233887
MT 0.9800422357402336
5-Block MinSLEM 0.8909885739345229
5-Block MinResist 0.8675200794119812
5-Block getMaxConnectivity 0.8694664535695392
5-Block Min |Z-W| 0.8841741753141091
3-Block MinSLEM 0.7141463212211673
3-Block MinResist 0.718263141982153
3-Block getMaxConnectivity 0.715183390367131
3-Block Min |Z-W| 0.7154940848780356
16
Full 0.4881613958305979
MT 0.982331559673419
8-Block MinSLEM 0.9617907227705366
8-Block MinResist 0.9432016873516004
8-Block getMaxConnectivity 0.9547829224218511
8-Block Min |Z-W| 0.984315427519989
4-Block MinSLEM 0.8305473874473303
4-Block MinResist 0.8108999202040502
4-Block getMaxConnectivity 0.8002305185657219
4-Block Min |Z-W| 0.9847811629632465
2-Block MinSLEM 0.5571562810349651
2-Block MinResist 0.42857145418996156
2-Block getMaxConnectivity 0.5061192608525976
2-Block Min |Z-W| 0.4285714564445944
17
Full 0.48901400167960424
MT 0.9842506483711779
18
Full 0.48976341389760036
MT 0.9858748663019703
9-Block MinSLEM 0.9703808981664126
9-Block MinResist 0.9543244228002966
9-Block getMaxConnectivity 0.9652620346221316
9-Block Min |Z-W| 0.9609559413299223
6-Block MinSLEM 0.9274863340746433
6-Block MinResist 0.904322100241284
6-Block getMaxConnectivity 0.9144368470036185
6-Block Min |Z-W| 0.9877641370733509
3-Block MinSLEM 0.7137591263145242
3-Block MinResist 0.7182631394317455
3-Block getMaxConnectivity 0.7156777756442414
3-Block Min |Z-W| 0.7151300210105594
2-Block MinSLEM 0.5571562751599648
2-Block MinResist 0.42857145889844994
2-Block getMaxConnectivity 0.5065757352882716
2-Block Min |Z-W| 0.428571459613919
19
Full 0.49042725685598576
MT 0.9872613997720258
20
Full 0.49101936604430657
MT 0.98845429318482
10-Block MinSLEM 0.9764918850332077
10-Block MinResist 0.9624236295054766
10-Block getMaxConnectivity 0.9726918301243737
10-Block Min |Z-W| 0.9899039312310299
5-Block MinSLEM 0.8909604820256735
5-Block MinResist 0.8675196430797856
5-Block getMaxConnectivity 0.8694657973902343
5-Block Min |Z-W| 0.884178006461683
4-Block MinSLEM 0.8306546228248091
4-Block MinResist 0.8108995981777251
4-Block getMaxConnectivity 0.8002313721274206
4-Block Min |Z-W| 0.9902135583151264
2-Block MinSLEM 0.5571562921308666
2-Block MinResist 0.42857146164342896
2-Block getMaxConnectivity 0.5073479496225984
2-Block Min |Z-W| 0.4285714623264254
21
Full 0.4915507533342952
MT 0.9894878242785429
7-Block MinSLEM 0.9487432477656033
7-Block MinResist 0.9275963069810571
7-Block getMaxConnectivity 0.9415195185314065
7-Block Min |Z-W| 0.9379767047811336
3-Block MinSLEM 0.7135615903082121
3-Block MinResist 0.7182631426601115
3-Block getMaxConnectivity 0.7154026655455449
3-Block Min |Z-W| 0.7151513652959529
22
Full 0.49203028591117437
MT 0.9903890671242299
11-Block MinSLEM 0.980965610799092
11-Block MinResist 0.9685141496628001
11-Block getMaxConnectivity 0.9778676717877701
11-Block Min |Z-W| 0.9734129802155151
2-Block MinSLEM 0.5571562777399365
2-Block MinResist 0.4285714662419582
2-Block getMaxConnectivity 0.5091397539683903
2-Block Min |Z-W| 0.428571466254424
23
Full 0.4924651929745758
MT 0.9911795563889043
24
Full 0.49286136387229174
MT 0.9918767125518874
12-Block MinSLEM 0.9842233582587118
12-Block MinResist 0.9732178088159736
12-Block getMaxConnectivity 0.9816646292439926
12-Block Min |Z-W| 0.9929855790344766
8-Block MinSLEM 0.9617899146861579
8-Block MinResist 0.9432016721407476
8-Block getMaxConnectivity 0.95478284148592
8-Block Min |Z-W| 0.9930940221197248
6-Block MinSLEM 0.9274864246060913
6-Block MinResist 0.904322406355109
6-Block getMaxConnectivity 0.9144367580545995
6-Block Min |Z-W| 0.9931276194555293
4-Block MinSLEM 0.8310241284105948
4-Block MinResist 0.8108982477563402
4-Block getMaxConnectivity 0.8002305221536402
4-Block Min |Z-W| 0.9931862824346592
3-Block MinSLEM 0.713933211150172
3-Block MinResist 0.7182631402120045
3-Block getMaxConnectivity 0.7152482065677889
3-Block Min |Z-W| 0.7152800293010432
2-Block MinSLEM 0.5571562808431788
2-Block MinResist 0.4285714310144653
2-Block getMaxConnectivity 0.5099386146815007
2-Block Min |Z-W| 0.42857147069660706

# Figure 9
df[df["obj"].isin(["MinResist", "Full", "MT"])].groupby(['obj', 'blocks'])['tau'].plot(marker='o', linestyle='-', figsize=(10, 6))
keys = ['Full', 'MT', '2-Block', '3-Block', '4-Block', '5-Block', '6-Block', '7-Block', '8-Block', '9-Block', '10-Block', '11-Block', '12-Block']
plt.legend(keys, bbox_to_anchor=(1, 1))

# Make x axis integer
plt.xticks(np.arange(4, 25, 2))
plt.xlabel(r"$n$")
plt.title("Min Resistance Contraction over $n$ Lipschitz Strongly Monotone Operators by Block Number")
if saving: plt.savefig("figs/fig_block_count.pdf", format="pdf", bbox_inches="tight")

# Figure 9 - Test function
def testMatricesProj(n, cases, names, lc=2, mu=1):
    ls = np.ones(n)*lc
    mus = np.ones(n)*mu
    ls[n-1] = np.inf
    mus[n-1] = 0
    taus = {}
    for i, (Z, M) in enumerate(cases):
        tau, _ = getReducedContractionOptGamma(Z, M, ls, mus)
        print(names[i], tau)
        taus[names[i]] = tau

    return taus

# Figure 9
projresults = []
for n in range(4,26): #t < 120:
    print(n)
    cases, names, objs, blocksizes, blocks = getTestMatricesBlocks(n)
    taus = testMatricesProj(n, cases, names)
    # Each key in taus is a column, the row is the value of n
    for i in range(len(names)):
        k = names[i]
        v = taus[k]
        
        projresults.append({"n": n, "name": k, "obj": objs[i], "blocksize": blocksizes[i], "blocks":blocks[i], "tau": v})

projdf = pd.DataFrame(projresults)
if saving: projdf.to_csv("figs/taus_proj.csv", index=False)
projdf.set_index('n', inplace=True)

4
Full 0.7521521274547781
MT 0.9177841369745806
2-Block MinSLEM 0.8191387832636592
2-Block MinResist 0.7870338922102957
2-Block getMaxConnectivity 0.8005884054570187
2-Block Min |Z-W| 0.7870338929300079
5
Full 0.804479212788819
MT 0.9475773165769151
6
Full 0.8385587406685645
MT 0.9636790407235943
3-Block MinSLEM 0.9045292083157187
3-Block MinResist 0.898926309334002
3-Block getMaxConnectivity 0.9031665202361303
3-Block Min |Z-W| 0.9058115649311607
2-Block MinSLEM 0.8708680920365539
2-Block MinResist 0.8616700832093386
2-Block getMaxConnectivity 0.8648381115685787
2-Block Min |Z-W| 0.861670081386794
7
Full 0.8625324131656764
MT 0.9733944272799407
8
Full 0.8803148080938081
MT 0.9796953548701692
4-Block MinSLEM 0.9477391187338774
4-Block MinResist 0.9366150047224902
4-Block getMaxConnectivity 0.9409099369418765
4-Block Min |Z-W| 0.9733653147481017
2-Block MinSLEM 0.8996656566051273
2-Block MinResist 0.8973584358184201
2-Block getMaxConnectivity 0.8970613583905448
2-Block Min |Z-W| 0.8973584363621916
9
Full 0.8940287042166508
MT 0.9840061423604585
3-Block MinSLEM 0.9310037419769321
3-Block MinResist 0.9226238048817228
3-Block getMaxConnectivity 0.9308683371516419
3-Block Min |Z-W| 0.9329275875236587
10
Full 0.9049261834616941
MT 0.987082059420084
5-Block MinSLEM 0.9690123236522284
5-Block MinResist 0.9580312775971841
5-Block getMaxConnectivity 0.962997391632803
5-Block Min |Z-W| 0.962987675833271
2-Block MinSLEM 0.9180267565103377
2-Block MinResist 0.9183673760071445
2-Block getMaxConnectivity 0.9166915390957012
2-Block Min |Z-W| 0.9183673754845927
11
Full 0.9137934450814728
MT 0.9893521326450776
12
Full 0.9211491280506132
MT 0.9910743537751785
6-Block MinSLEM 0.9794380910215909
6-Block MinResist 0.9701555310272644
6-Block getMaxConnectivity 0.9756973867550393
6-Block Min |Z-W| 0.99010537100968
4-Block MinSLEM 0.9594824501758149
4-Block MinResist 0.9489711986483735
4-Block getMaxConnectivity 0.9564358341424667
4-Block Min |Z-W| 0.9875641253599383
3-Block MinSLEM 0.9461116981534052
3-Block MinResist 0.9376136641003483
3-Block getMaxConnectivity 0.9476885966471409
3-Block Min |Z-W| 0.9478983876056438
2-Block MinSLEM 0.9307344399513203
2-Block MinResist 0.9322246586548163
2-Block getMaxConnectivity 0.9300141956611677
2-Block Min |Z-W| 0.932224658678786
13
Full 0.9273491859704465
MT 0.992411364833255
14
Full 0.932645842257762
MT 0.993469781964156
7-Block MinSLEM 0.9854895007855449
7-Block MinResist 0.9778942634910617
7-Block getMaxConnectivity 0.9826255537601088
7-Block Min |Z-W| 0.9813969384739117
2-Block MinSLEM 0.9400431051319242
2-Block MinResist 0.942055224496676
2-Block getMaxConnectivity 0.9396299797588917
2-Block Min |Z-W| 0.9420552246279732
15
Full 0.937223186841485
MT 0.9943217922540158
5-Block MinSLEM 0.974590401353922
5-Block MinResist 0.9649800729630257
5-Block getMaxConnectivity 0.9710959252251317
5-Block Min |Z-W| 0.9698525014507775
3-Block MinSLEM 0.9561864907177353
3-Block MinResist 0.9478753095789055
3-Block getMaxConnectivity 0.9572607197744706
3-Block Min |Z-W| 0.9575677751989082
16
Full 0.9412182809698121
MT 0.9950176830138516
8-Block MinSLEM 0.989233292652846
8-Block MinResist 0.9829428678589387
8-Block getMaxConnectivity 0.9871235601177183
8-Block Min |Z-W| 0.9951900480677953
4-Block MinSLEM 0.9670904965450355
4-Block MinResist 0.9575260060505665
4-Block getMaxConnectivity 0.9660447128654579
4-Block Min |Z-W| 0.992835672581666
2-Block MinSLEM 0.9471521763588353
2-Block MinResist 0.9493931058277882
2-Block getMaxConnectivity 0.9469224619297711
2-Block Min |Z-W| 0.9493931076789691
17
Full 0.9447354996217051
MT 0.9955932803549359
18
Full 0.947855785340868
MT 0.996074790975494
9-Block MinSLEM 0.9917523326892693
9-Block MinResist 0.9865002011040845
9-Block getMaxConnectivity 0.9901485632323862
9-Block Min |Z-W| 0.9891810853182708
6-Block MinSLEM 0.9825422672544786
6-Block MinResist 0.9744372782097804
6-Block getMaxConnectivity 0.9802334300928883
6-Block Min |Z-W| 0.9955573342110063
3-Block MinSLEM 0.9625973759072839
3-Block MinResist 0.9552965512070167
3-Block getMaxConnectivity 0.9643509557610295
3-Block Min |Z-W| 0.9639152802684406
2-Block MinSLEM 0.9527572967367829
2-Block MinResist 0.955080203083856
2-Block getMaxConnectivity 0.9526415889171332
2-Block Min |Z-W| 0.9550802035807573
19
Full 0.9506426888675692
MT 0.996481611885768
20
Full 0.9531469126958941
MT 0.9968284121901984
10-Block MinSLEM 0.9933658104372783
10-Block MinResist 0.989021596450872
10-Block getMaxConnectivity 0.9921891728606781
10-Block Min |Z-W| 0.9972491010019479
5-Block MinSLEM 0.9784205059914235
5-Block MinResist 0.9700027732901295
5-Block getMaxConnectivity 0.9767083161681916
5-Block Min |Z-W| 0.9748735206170638
4-Block MinSLEM 0.9721919454461948
4-Block MinResist 0.9637516383754409
4-Block getMaxConnectivity 0.9724264674719896
4-Block Min |Z-W| 0.995348815017694
2-Block MinSLEM 0.9572893031112173
2-Block MinResist 0.9596176455086766
2-Block getMaxConnectivity 0.9572470127241994
2-Block Min |Z-W| 0.9596176491745421
21
Full 0.9554093767969243
MT 0.9971264271830746
7-Block MinSLEM 0.9874743330706199
7-Block MinResist 0.9805966970219246
7-Block getMaxConnectivity 0.9859754616666889
7-Block Min |Z-W| 0.983739405930636
3-Block MinSLEM 0.9674384858406184
3-Block MinResist 0.9608952673738557
3-Block getMaxConnectivity 0.9689905900505172
3-Block Min |Z-W| 0.9688040622222985
22
Full 0.9574634693504498
MT 0.9973843798894533
11-Block MinSLEM 0.9945861578488714
11-Block MinResist 0.9908850425394582
11-Block getMaxConnectivity 0.9936949377169232
11-Block Min |Z-W| 0.9930500080406014
2-Block MinSLEM 0.9610289962695155
2-Block MinResist 0.9633222012558504
2-Block getMaxConnectivity 0.9610325830839557
2-Block Min |Z-W| 0.9633222012846975
23
Full 0.9593367015363286
MT 0.997609135727751
24
Full 0.9610519693853471
MT 0.9978061502071545
12-Block MinSLEM 0.9954966213190756
12-Block MinResist 0.992309516621772
12-Block getMaxConnectivity 0.9947939665976838
12-Block Min |Z-W| 0.9982502783079058
8-Block MinSLEM 0.9904306964795818
8-Block MinResist 0.9847609731860296
8-Block getMaxConnectivity 0.9894469413820468
8-Block Min |Z-W| 0.9978752519048323
6-Block MinSLEM 0.9855246976114859
6-Block MinResist 0.9776570058927866
6-Block getMaxConnectivity 0.9841645839948792
6-Block Min |Z-W| 0.9974805118666883
4-Block MinSLEM 0.9756430722135762
4-Block MinResist 0.9684466642076521
4-Block getMaxConnectivity 0.9763035559878469
4-Block Min |Z-W| 0.996738808067415
3-Block MinSLEM 0.9715913532679965
3-Block MinResist 0.9652615861029722
3-Block getMaxConnectivity 0.972578833986059
3-Block Min |Z-W| 0.9726116359656548
2-Block MinSLEM 0.9641672314385882
2-Block MinResist 0.96640390378116
2-Block getMaxConnectivity 0.9642023236849637
2-Block Min |Z-W| 0.9664039049574353
25
Full 0.9626284123300337
MT 0.9979798028060392
5-Block MinSLEM 0.9811923650459194
5-Block MinResist 0.9738455481222618
5-Block getMaxConnectivity 0.9805674471032336
5-Block Min |Z-W| 0.9784107033246009

# Figure 9 Class 2 

projdf[projdf["obj"].isin(["MinResist", "Full", "MT"])].groupby(['obj', 'blocks'])['tau'].plot(marker='o', linestyle='-', figsize=(10, 6))

keys = ['Full', 'MT', '2-Block', '3-Block', '4-Block', '5-Block', '6-Block', '7-Block', '8-Block', '9-Block', '10-Block', '11-Block', '12-Block'] # If using df directly
plt.legend(keys, bbox_to_anchor=(1, 1))
plt.title("Min Resistance Contraction over $n-1$ Lipschitz Strongly Monotone Operators")
if saving: plt.savefig("figs/fig_block_count_proj.pdf", format="pdf", bbox_inches="tight")

Figure 10 combines the iteration time and contraction factor to depict how removing edges impacts algorithm design. We see that the 2-Block design strikes the optimal balance between connectivity (for faster contraction) and sparseness (for parallelism and faster iteration times)

This uses the getIterationTime function to return the iteration time.

# Figure 10 - pareto curves

from matplotlib.colors import CenteredNorm
def getTaus(Z, W, l=None, mu=None):

    M = getMfromWCholesky(W)
    h = getReducedContractionFactor(Z, M, l, mu, gamma=0.5)
    m, _ = getReducedContractionOptGamma(Z, M, l, mu)
    return h, m

def getContractTime(n, l, mu, builder=None, tol=0.5):
    Z_fixed = {}

    t = np.ones(n)
    ll = np.ones((n,n))
    results = {}
    for i in range(n//2):
        print(i)
        start = 0
        for k in range(2):
            for j in range(start, i + k*(n//2)):
                Z_fixed[(i + k*(n//2),j)] = 0
            Zs, Ws = builder(n, fixed_Z=Z_fixed)
            h, m = getTaus(Zs, Ws, l, mu)
            tt, _, _ = getIterationTime(t, ll, Zs, Ws, itrs=2)
            # r is number of entries in Z_fixed
            r = len(Z_fixed)
            tif = np.ceil(np.log(tol)/np.log(h)) # Total iterations fixed
            tio = np.ceil(np.log(tol)/np.log(m)) # Total iterations opt step
            ttf = tif*tt
            tto = tio*tt
            results[r] = {'Fixed':h , 'M Opt': m, 'Time': tt, 'Removed Edges': r, 'Fixed Total Iterations': tif, 'Opt Step Total Iterations': tio, 'Fixed Total Time': ttf, 'Opt Step Total Time': tto}
            print(results[r])
            start = n//2

    counter = 0
    for i in range(n//2 - 2):
        for j in range(i+1):
            r = n//2 + j
            c = n//2 - 1 - i + j
            Z_fixed[(r, c)] = 0
            counter += 1
            if counter % 10 == 0:
                Ls, Ws = builder(n, fixed_Z=Z_fixed)
                h, m = getTaus(Ls, Ws, l, mu)
                tt = 4
                # r is number of entries in Z_fixed
                rem = len(Z_fixed)    
                tif = np.ceil(np.log(tol)/np.log(h)) # Total iterations fixed
                tio = np.ceil(np.log(tol)/np.log(m)) # Total iterations opt step
                ttf = tif*tt
                tto = tio*tt
                results[rem] = {'Fixed':h , 'M Opt': m, 'Time': tt, 'Removed Edges': rem, 'Fixed Total Iterations': tif, 'Opt Step Total Iterations': tio, 'Fixed Total Time': ttf, 'Opt Step Total Time': tto}
                print(results[rem])
    # Create dataframe of taus and times
    df = pd.DataFrame(results)

    return df

# Figure 10
n = 24
l = np.ones(n)*2
l[n-1] = np.inf
mu = np.ones(n)
mu[n-1]= 0

df = getContractTime(n, l, mu, builder=getMaxConnectivity)
df2 = df.T

0
{'Fixed': 0.9798415277460015, 'M Opt': 0.9610519693852897, 'Time': 47.99999999999555, 'Removed Edges': 0, 'Fixed Total Iterations': 35.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1679.9999999998442, 'Opt Step Total Time': 863.99999999992}
{'Fixed': 0.9798415277460015, 'M Opt': 0.9610519693852897, 'Time': 47.99999999999555, 'Removed Edges': 0, 'Fixed Total Iterations': 35.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1679.9999999998442, 'Opt Step Total Time': 863.99999999992}
1
{'Fixed': 0.9802934217043753, 'M Opt': 0.9611344874981325, 'Time': 46.00000000002653, 'Removed Edges': 1, 'Fixed Total Iterations': 35.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1610.0000000009286, 'Opt Step Total Time': 828.0000000004776}
{'Fixed': 0.9802906549740917, 'M Opt': 0.9612008911742602, 'Time': 44.00000000001397, 'Removed Edges': 2, 'Fixed Total Iterations': 35.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1540.0000000004889, 'Opt Step Total Time': 792.0000000002515}
2
{'Fixed': 0.9803115291134975, 'M Opt': 0.9612841961148996, 'Time': 41.999999999993534, 'Removed Edges': 4, 'Fixed Total Iterations': 35.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1469.9999999997738, 'Opt Step Total Time': 755.9999999998836}
{'Fixed': 0.9803217937756855, 'M Opt': 0.9613404323358163, 'Time': 40.00000000033529, 'Removed Edges': 6, 'Fixed Total Iterations': 35.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1400.0000000117352, 'Opt Step Total Time': 720.0000000060352}
3
{'Fixed': 0.9803594774713668, 'M Opt': 0.9614511691184635, 'Time': 38.00000000019668, 'Removed Edges': 9, 'Fixed Total Iterations': 35.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1330.0000000068837, 'Opt Step Total Time': 684.0000000035402}
{'Fixed': 0.9803725407513265, 'M Opt': 0.9615160721244879, 'Time': 36.00000000005643, 'Removed Edges': 12, 'Fixed Total Iterations': 35.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1260.0000000019752, 'Opt Step Total Time': 648.0000000010158}
4
{'Fixed': 0.9804160035586761, 'M Opt': 0.9616398317573042, 'Time': 34.00000000000089, 'Removed Edges': 16, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1224.000000000032, 'Opt Step Total Time': 612.000000000016}
{'Fixed': 0.9803159903286937, 'M Opt': 0.9616501789931347, 'Time': 32.00000000018983, 'Removed Edges': 20, 'Fixed Total Iterations': 35.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1120.000000006644, 'Opt Step Total Time': 576.000000003417}
5
{'Fixed': 0.9805064778081075, 'M Opt': 0.9618143984054425, 'Time': 29.99999999993441, 'Removed Edges': 25, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1079.9999999976387, 'Opt Step Total Time': 539.9999999988194}
{'Fixed': 0.9805369849103691, 'M Opt': 0.9618911233999475, 'Time': 28.000000000015845, 'Removed Edges': 30, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 1008.0000000005705, 'Opt Step Total Time': 504.00000000028524}
6
{'Fixed': 0.980558744118224, 'M Opt': 0.9619696358581257, 'Time': 26.00000000018838, 'Removed Edges': 36, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 936.0000000067816, 'Opt Step Total Time': 468.0000000033908}
{'Fixed': 0.9806165695795973, 'M Opt': 0.9620812351429494, 'Time': 23.999999999924594, 'Removed Edges': 42, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 18.0, 'Fixed Total Time': 863.9999999972854, 'Opt Step Total Time': 431.9999999986427}
7
{'Fixed': 0.9805803035543116, 'M Opt': 0.9622244218633337, 'Time': 22.000000000001247, 'Removed Edges': 49, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 19.0, 'Fixed Total Time': 792.0000000000449, 'Opt Step Total Time': 418.0000000000237}
{'Fixed': 0.980609771943829, 'M Opt': 0.9623042101456656, 'Time': 19.999999999880224, 'Removed Edges': 56, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 19.0, 'Fixed Total Time': 719.9999999956881, 'Opt Step Total Time': 379.9999999977243}
8
{'Fixed': 0.9806705602427209, 'M Opt': 0.962424427613141, 'Time': 17.99999999994352, 'Removed Edges': 64, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 19.0, 'Fixed Total Time': 647.9999999979667, 'Opt Step Total Time': 341.99999999892685}
{'Fixed': 0.9807014558127024, 'M Opt': 0.9624855242813694, 'Time': 15.999999999994683, 'Removed Edges': 72, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 19.0, 'Fixed Total Time': 575.9999999998086, 'Opt Step Total Time': 303.999999999899}
9
{'Fixed': 0.9808093208410873, 'M Opt': 0.9627105025538466, 'Time': 14.000000000000362, 'Removed Edges': 81, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 19.0, 'Fixed Total Time': 504.0000000000131, 'Opt Step Total Time': 266.0000000000069}
{'Fixed': 0.9808410857473983, 'M Opt': 0.9627679103019492, 'Time': 11.9999999999999, 'Removed Edges': 90, 'Fixed Total Iterations': 36.0, 'Opt Step Total Iterations': 19.0, 'Fixed Total Time': 431.9999999999964, 'Opt Step Total Time': 227.99999999999812}
10
{'Fixed': 0.9810231148410529, 'M Opt': 0.96304306184268, 'Time': 10.000000000165414, 'Removed Edges': 100, 'Fixed Total Iterations': 37.0, 'Opt Step Total Iterations': 19.0, 'Fixed Total Time': 370.00000000612033, 'Opt Step Total Time': 190.00000000314287}
{'Fixed': 0.9810393061615912, 'M Opt': 0.9630413776592729, 'Time': 7.999999999999933, 'Removed Edges': 110, 'Fixed Total Iterations': 37.0, 'Opt Step Total Iterations': 19.0, 'Fixed Total Time': 295.99999999999756, 'Opt Step Total Time': 151.99999999999872}
11
{'Fixed': 0.9813641777345405, 'M Opt': 0.9642038633930519, 'Time': 3.9999999999999907, 'Removed Edges': 121, 'Fixed Total Iterations': 37.0, 'Opt Step Total Iterations': 20.0, 'Fixed Total Time': 147.99999999999966, 'Opt Step Total Time': 79.99999999999982}
{'Fixed': 0.9813651787490266, 'M Opt': 0.9642023236823605, 'Time': 3.9999999999999907, 'Removed Edges': 132, 'Fixed Total Iterations': 37.0, 'Opt Step Total Iterations': 20.0, 'Fixed Total Time': 147.99999999999966, 'Opt Step Total Time': 79.99999999999982}
{'Fixed': 0.9838400213581495, 'M Opt': 0.9652843225004955, 'Time': 4, 'Removed Edges': 142, 'Fixed Total Iterations': 43.0, 'Opt Step Total Iterations': 20.0, 'Fixed Total Time': 172.0, 'Opt Step Total Time': 80.0}
{'Fixed': 0.9853883421029684, 'M Opt': 0.9669002363662484, 'Time': 4, 'Removed Edges': 152, 'Fixed Total Iterations': 48.0, 'Opt Step Total Iterations': 21.0, 'Fixed Total Time': 192.0, 'Opt Step Total Time': 84.0}
{'Fixed': 0.9868013528993385, 'M Opt': 0.968647136449619, 'Time': 4, 'Removed Edges': 162, 'Fixed Total Iterations': 53.0, 'Opt Step Total Iterations': 22.0, 'Fixed Total Time': 212.0, 'Opt Step Total Time': 88.0}
{'Fixed': 0.9886689502714882, 'M Opt': 0.9699045194765014, 'Time': 4, 'Removed Edges': 172, 'Fixed Total Iterations': 61.0, 'Opt Step Total Iterations': 23.0, 'Fixed Total Time': 244.0, 'Opt Step Total Time': 92.0}
{'Fixed': 0.9912682146833666, 'M Opt': 0.9724128400046071, 'Time': 4, 'Removed Edges': 182, 'Fixed Total Iterations': 80.0, 'Opt Step Total Iterations': 25.0, 'Fixed Total Time': 320.0, 'Opt Step Total Time': 100.0}

if saving: df.to_csv('figs/n24.csv', index=False)

# Plot total time vs removed edges
plt.plot(df.iloc[3], df.iloc[6], df.iloc[3], df.iloc[7])
# Subtitle
plt.suptitle('Time Units to 0.5 Contraction vs Removed Edges, $n=$'+str(n))


plt.legend( ['Fixed Step Size', 'Optimal Step Size'])
plt.xlabel('Removed Edges')
plt.ylabel('Time Units')
if saving: plt.savefig("figs/fig_optimal_sparsity.pdf", format="pdf", bbox_inches="tight")
plt.show()

# Figure 10
# Plot the top three rows of df as a function of time with points annotated
plt.plot(df.iloc[2], df.iloc[0], marker='<')
plt.plot(df.iloc[2], df.iloc[1], marker='<')

plt.title('Contraction vs Time as Edges are Removed')
# Add legend
plt.legend( ['Fixed Step Size', 'Optimal Step Size'])
# Annotate points with number of removed edges
lx, ly = (0,0)
right = 0
up = 10
for j in range(len(df.columns)):
        print(j,(df.iloc[2,j], df.iloc[1,j]))
        if df.iloc[2,j] < 4:
                right = 10
                up = 1
        if df.iloc[2,j] == 4.0:
                up = 0
        if abs(lx - df.iloc[2,j]) > 1 or abs(ly - df.iloc[1,j]) > 0.001:
                plt.annotate(df.columns[j], (df.iloc[2,j], df.iloc[0,j]), textcoords="offset points", xytext=(right,up), ha='center', fontsize=8)
                plt.annotate(df.columns[j], (df.iloc[2,j], df.iloc[1,j]), textcoords="offset points", xytext=(right,up), ha='center', fontsize=8)
                lx = df.iloc[2,j]
                ly = df.iloc[1,j]
        
# Annotate with points

plt.xlabel('Average Iteration Time')
plt.ylabel(r'Contraction Rate ($\tau$)')

if saving: plt.savefig("figs/fig_pareto.pdf", format="pdf", bbox_inches="tight")
plt.show()

0 (47.99999999999555, 0.9610519693852897)
1 (46.00000000002653, 0.9611344874981325)
2 (44.00000000001397, 0.9612008911742602)
3 (41.999999999993534, 0.9612841961148996)
4 (40.00000000033529, 0.9613404323358163)
5 (38.00000000019668, 0.9614511691184635)
6 (36.00000000005643, 0.9615160721244879)
7 (34.00000000000089, 0.9616398317573042)
8 (32.00000000018983, 0.9616501789931347)
9 (29.99999999993441, 0.9618143984054425)
10 (28.000000000015845, 0.9618911233999475)
11 (26.00000000018838, 0.9619696358581257)
12 (23.999999999924594, 0.9620812351429494)
13 (22.000000000001247, 0.9622244218633337)
14 (19.999999999880224, 0.9623042101456656)
15 (17.99999999994352, 0.962424427613141)
16 (15.999999999994683, 0.9624855242813694)
17 (14.000000000000362, 0.9627105025538466)
18 (11.9999999999999, 0.9627679103019492)
19 (10.000000000165414, 0.96304306184268)
20 (7.999999999999933, 0.9630413776592729)
21 (3.9999999999999907, 0.9642038633930519)
22 (3.9999999999999907, 0.9642023236823605)
23 (4.0, 0.9652843225004955)
24 (4.0, 0.9669002363662484)
25 (4.0, 0.968647136449619)
26 (4.0, 0.9699045194765014)
27 (4.0, 0.9724128400046071)

Figure 11 just shows how we removed the edges

# Figure 11 - edge deletion ordering

def cplot(M, c, title=''):
    vmin = 1
    vmax = np.max(ordermat)
    cm = plt.cm.viridis #coolwarm
    fig, ax = plt.subplots(1, 1, figsize=(12, 6))
    x = np.arange(n)
    y = -np.arange(n)
    pc = ax.pcolormesh(x, y, M, cmap=cm , vmin=vmin, vmax=vmax)
    ax.set_title(title) #'Edge Removal Ordering'
    ax.axis('off')
    ax.set_aspect('equal')

    # Add colorbar
    plt.tight_layout()
    
    cbar = fig.colorbar(pc) 
    cbar.set_ticks(np.arange(1, c, 4)) 
    return fig

# Display edge removal ordering
def order_graph(n):
    order_matrix = np.zeros((n,n)) 
    counter = 1
    for i in range(n//2):
        start = 0
        for k in range(2):
            for j in range(start, i + k*(n//2)):
                r = i + k*(n//2)
                c = j
                order_matrix[(r,j)] = counter
                order_matrix[(c,r)] = counter
            
            counter += 1
            start = n//2

    ccounter = 0
    for i in range(n//2 - 2):
        for j in range(i+1):
            r = n//2 + j
            c = n//2 - 1 - i + j
            order_matrix[(r, c)] = counter
            order_matrix[(c, r)] = counter
            ccounter += 1
            if ccounter % 10 == 0:
                counter += 1
            
    return order_matrix, counter
ordermat, c = order_graph(24)
ordermat = np.ma.masked_array(ordermat, ordermat < 1)

fig = cplot(ordermat, c, 'Edge Removal Ordering')
if saving: fig.savefig("figs/fig_edge_removal.pdf")

Figure 12 then combines the iteration time and contraction factor to find the total time required for a given contraction level, and compares the 2-Block design with the constrained minimum iteration time design, using the minimum total effective resistance objective for each.

It demonstrates some additional parameters in the getMinIteration function. minW=n//2 sets the minimum number of edges for each node in the graph of $W$ to be $\frac{n}{2}$, which we call the constrained minimum iteration time design.

We also introduce the getMetrics function, which returns the iteration time, contraction factor, total cycles, and total time for a given contraction target, over a given design ($Z, W$), compute and communication times ($t, l$), and operator parameters (using the default values in Class 1).

# Figure 12
block_results = []
min_results = []
for n in range(6,10):
    for i in range(40):
        t = np.random.rand(n)*10+1

        # Communication time
        l = np.random.rand(n, n)*1.5 + 0.5
        l = np.tril(l, -1) + np.tril(l, -1).T
        Zb, Wb = getBlockMin(n, n//2, builder=getMinResist)
        c, tau, cc, tt = getMetrics(Zb, Wb, t, l,contraction_target=0.01)
        block_results.append((n, i, c, tau, cc, tt))
        Z, W = getMinIteration(n, builder=getMinResist, t=t, l=l, c=0.01, r=n, minW=n//2)
        if Z is not None:
            c, tau, cc, ttc = getMetrics(Z, W, t, l,contraction_target=0.01)
        else:
            c, tau, cc, ttc = (None, None, None, None)
        min_results.append((n, i, c, tau, cc, ttc))
        print(n,i,tt,ttc)
        

6 0 670.4886897131029 670.4886897131029
6 1 778.643903904538 694.5692581525664
6 2 615.006840052717 615.006840052717
6 3 664.7031827966202 664.7031827966202
6 4 698.4652093107794 609.2499864585676
6 5 688.0880207701939 688.0880207701939
6 6 577.4186832479074 513.2587602400457
6 7 675.5954285310031 675.5954285310031
6 8 652.2238505242201 632.6068381216719
6 9 623.2432588359057 498.47284290251457
6 10 530.4714698292019 530.4714698292019
6 11 487.89630928740723 487.89630928740723
6 12 616.5828737639896 616.5828737639896
6 13 545.1163811472464 496.95239405209304
6 14 715.8385535884668 544.7866595980308
6 15 718.7565870708886 718.7565870708886
6 16 508.98055829220084 508.98055829220084
6 17 618.5273437625365 618.5273437625365
6 18 587.1497931809487 587.1497931809487
6 19 697.2152163421836 697.2152163421836
6 20 525.217813021972 525.217813021972
6 21 653.3343949915621 653.3343949915621
6 22 451.8401413427987 395.53697807569
6 23 541.1292509237803 541.1292509237803
6 24 560.1621320320938 560.1621320320938
6 25 627.8152551777415 513.6085171905685
6 26 697.0303838769128 415.72627354570216
6 27 591.7380582303297 591.7380582303297
6 28 643.5697463454466 579.1653239472738
6 29 570.225004093676 570.225004093676
6 30 670.7815318254211 553.1210091819509
6 31 651.4000885000901 651.4000885000901
6 32 590.7057065157435 590.7057065157435
6 33 501.7307153879898 501.7307153879898
6 34 562.4305705164893 562.4305705164893
6 35 610.6595096137099 610.6595096137099
6 36 571.2032212576437 464.3189390758997
6 37 517.0737763095983 466.544290865819
6 38 377.0452637567923 351.44461502719867
6 39 689.041606842278 577.3529920169221
7 0 981.2245958002236 786.427289915642
7 1 940.6699817649405 647.5730411744605
7 2 925.0920189003568 785.2373220527802
7 3 857.3204311391953 745.4147206690646
7 4 982.8055531664709 870.8995244936355
7 5 1053.766136245001 724.8240692802201
7 6 1129.0206508574934 769.8780816183216
7 7 1091.3057839459616 1024.4911440881733
7 8 626.6460464052768 562.7025723831766
7 9 997.8021089014212 685.620755299691
7 10 797.2032611075188 585.7003551257465
7 11 885.0031300908411 848.880553339001
7 12 1058.9747887818658 706.7751498158004
7 13 842.7418733900938 619.1572947659758
7 14 818.8331804928644 718.5678931618287
7 15 987.3290374765884 866.4316043159722
7 16 1064.5623716096825 719.9327341961316
7 17 976.8999143479339 671.1810208221348
7 18 914.6632549494241 1119.9958223882777
7 19 808.3796553566091 530.2767466364537
7 20 1034.275498773276 759.8758766500566
7 21 927.6688513933602 692.8621625075766
7 22 1060.2882154509803 784.9831485046384
7 23 1012.2134651303912 777.0082395230768
7 24 1195.958656756382 1003.1196132371452
7 25 1148.8716751623185 645.0757853899174
7 26 1115.9148391788121 670.7122101748075
7 27 1029.6629370800192 792.2840369171759
7 28 880.0399189098284 622.0091332190075
7 29 907.056906275165 795.988713669901
7 30 1096.0949150977958 777.5918321005179
7 31 840.4031523977338 617.4390507417304
7 32 920.9986365840482 808.2232933097077
7 33 1056.4003240787945 1293.5514171010198
7 34 1147.7506188719708 750.0765387475235
7 35 1075.8273218346242 844.2669768069983
7 36 955.2703771781606 734.3842685437616
7 37 1115.1376903156365 577.3904988889897
7 38 980.5479232805474 694.7159873902991
7 39 1131.1777553101992 725.4714377472677
8 0 911.6232630181286 810.4387510039699
8 1 891.9799263850759 891.9799263850759
8 2 920.179239297414 920.179239297414
8 3 800.2383349299959 800.2383349299959
8 4 753.3587758016674 753.3587758016674
8 5 845.0517049325248 845.0517049325248
8 6 912.1179851660208 912.1179851660208
8 7 848.446268018055 848.446268018055
8 8 982.3775931007713 982.3775931007713
8 9 952.1745880593768 952.1745880593768
8 10 950.707346055344 950.707346055344
8 11 829.9867654940377 633.3932317934864
8 12 969.4493413955097 736.1760875191544
8 13 997.2671015244088 864.122605905095
8 14 1018.3202726756092 1018.3202726756092
8 15 1046.8438441968804 1046.8438441968804
8 16 861.2354279979197 861.2354279979197
8 17 811.3450272417898 811.3450272417898
8 18 1036.2270906794056 821.6057232227172
8 19 939.266566851289 835.6024580003277
8 20 994.7175797539742 765.0848878706779
8 21 877.9064353334034 877.9064353334034
8 22 982.6989485839769 982.6989485839769
8 23 991.3768071185148 755.2427074708362
8 24 855.6159120155098 855.6159120155098
8 25 982.8271720948637 882.7642218776714
8 26 829.3763935963367 829.3763935963367
8 27 893.5883959113827 676.0339410506026
8 28 865.4880113161685 865.4880113161685
8 29 892.4905049124536 794.4842759862287
8 30 971.8327334063216 800.8026999372078
8 31 1008.9261947189451 826.3379529075881
8 32 894.7718786666485 752.2926906433577
8 33 815.7798898160743 815.7798898160743
8 34 992.1038858840896 992.1038858840896
8 35 691.1516078532711 691.1516078532711
8 36 947.354286830722 947.354286830722
8 37 922.6948858772778 922.6948858772778
8 38 918.150100608076 728.1685261482883
8 39 947.8385952610337 774.5080911243963
9 0 1478.2742502043352 979.6073962049029
9 1 1093.6531782621405 684.9243844953108
9 2 1201.525123664246 828.9801754145494
9 3 1364.5336587289812 950.0396165707325
9 4 1413.3915046389948 1005.4663986748551
9 5 974.688508500644 640.1167130385393
9 6 1350.6261230564066 862.4545701796865
9 7 1162.8223692152496 1016.324466625392
9 8 1210.9060082553626 806.1813811375092
9 9 1423.4552052809431 1103.312451259093
9 10 1291.5338676719557 791.3902376658899
9 11 1205.7916929464434 975.5444833475439
9 12 1451.529222078253 1107.8159920741818
9 13 1386.8755933773932 1005.4398941073596
9 14 1258.419176099138 912.2887549393351
9 15 1392.3203182854224 928.9888255176635
9 16 1360.4958752153861 1070.5541313171288
9 17 1037.1756829668498 793.579455126998
9 18 1294.3366325644306 1121.158048648846
9 19 1311.584916674905 987.2253639656058
9 20 1203.5323775901584 898.7089991280843
9 21 1129.3903829021258 1203.448768732348
9 22 1307.7908372303818 1082.822479286563
9 23 1271.0210381687546 796.4761526137635
9 24 1327.68395677353 759.956492346676
9 25 1290.9681193745478 836.009573751264
9 26 1312.1585307875869 815.8914406834145
9 27 1017.3591417545098 574.882217494349
9 28 912.1199881927741 801.5857158873254
9 29 1350.7925484160241 969.13105880696
9 30 1402.3920226256741 1075.8709172135295
9 31 1247.046436299857 975.9273880798273
9 32 1236.073399741868 720.2962432598588
9 33 1413.9913547720375 1159.009308067732
9 34 1430.4663131494335 1004.469912472709
9 35 1368.3719247419199 1054.319352635388
9 36 1149.2680010915517 906.3371960222526
9 37 1293.6039989563678 938.7706549615914
9 38 1402.9186006253822 1075.4741357450423
9 39 1144.0111659936217 849.204131805314

# Figure 12
block_df = pd.DataFrame(block_results)
min_df = pd.DataFrame(min_results)

if saving: 
    block_df.to_csv('figs/block_df2.csv')
    min_df.to_csv('figs/min_cycle_min_edge.csv')

# Figure 12
# Paired density chart for block vs min cycle

fig, axs = plt.subplots(2, 2, sharey=True)
n = 6
key = 5 #'5' #use string if loading from csv
start = 0
for i in range(2):
    for j in range(2):
        end = 40 + start
        axs[i,j] = block_df[key][start:end].plot.kde(ax=axs[i,j], color='blue')
        min_df[key][start:end].plot.kde(ax=axs[i,j], color='red')
        start += 40
        axs[i,j].set_title(r"$n=$"+str(n))
        n += 1

fig.legend([r'$d$-Block Design', 'Con. Min Iter Time'], loc=(.65,.8))

for ax in axs.flat:
    ax.set(xlabel='Time')
for ax in fig.get_axes():
    ax.label_outer()
if saving: fig.savefig("figs/fig_constrained_mincycle_histogram.pdf", format="pdf", bbox_inches="tight")