Network Optimization Algorithms

---

# 🌊 Network Optimization Algorithms

### C++ Implementations in `digraphx-cpp` & `netoptim-cpp`

@luk036 👨‍💻 · 2026 📅

---

### 📋 Agenda

#### 1️⃣ Foundation
- What is Network Optimization?
- Negative Cycle Detection (Howard's Algorithm)

#### 2️⃣ Core Algorithms
- Minimum Cycle Ratio (MCR)
- Parametric Network Problems

]

#### 3️⃣ Cutting-Plane Methods
- NetworkOracle Architecture
- Separation Oracles

#### 4️⃣ Applications
- Optimal Matrix Scaling ✨
- Clock Skew Scheduling ⏰
- Performance Analysis 📊

]

---

## ⚙️ Foundation: Negative Cycles

---

### 🔍 What is a Negative Cycle? 🔄

A **negative cycle** is a directed cycle where the sum of edge weights is negative:

$$ \sum_{(i,j) \in {\color{lime}C}} {\color{coral}w}(i,j) < 0 $$

<pre>
graph LR
    A[("v₁")] -->|"-2"| B[("v₂")]
    B -->|"1"| C[("v₃")]
    C -->|"1"| A
    style A fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style B fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style C fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
</pre>

]

This 3-cycle has total weight $-2 + 1 + 1 = 0$ (non-negative).

---

### Why Negative Cycles Matter ❓

In network optimization, negative cycles indicate **infeasibility**:

- **Shortest paths**: A negative cycle means no finite shortest path exists
- **Feasibility constraints**: ${\color{firebrick}u_j} - {\color{firebrick}u_i} \le {\color{coral}w}(i,j)$ — a negative cycle means the system has no solution
- **Cutting-plane methods**: Finding negative cycles generates **separation cuts**

<pre>
graph LR
    subgraph Input["Input Constraints"]
        A1("Edge Weights\nw(i,j)")
    end
    subgraph Detection["Cycle Detection"]
        B1{"Negative\nCycle?"}
    end
    subgraph Output["Result"]
        C1("Feasible ✅")
        C2("Infeasible ❌\n→ Cutting Plane")
    end
    Input --> Detection
    Detection -->|"No"| C1
    Detection -->|"Yes"| C2
    style Input fill:#e3f2fd,stroke:#1565c0,stroke-width:3px
    style Detection fill:#fff3e0,stroke:#e65100,stroke-width:3px
    style Output fill:#e8f5e9,stroke:#2e7d32,stroke-width:3px
</pre>

]

---

### Howard's Algorithm ⏱️

Howard's **policy iteration** algorithm for negative cycle detection:

1. **Relax**: Perform Bellman-Ford style edge relaxations
2. **Detect**: Find cycles in the predecessor graph
3. **Verify**: Confirm the cycle is actually negative
4. **Yield**: Return found cycles via generator

.font-sm.mb-xs[
```cpp
NegCycleFinder finder(graph);

// Generator-based enumeration
for (const auto& cycle : finder.howard(dist, get_weight)) {
    // Process the negative cycle
}
```
]

**Why Howard's?** 🏆

- No source node required
- Faster than Bellman-Ford for dense graphs
- **Generator-based** — memory-efficient cycle enumeration
- Reuses distance estimates between iterations

---

### Policy Iteration Visualized 🎯

<pre>
graph LR
    Start("Initialize\ndist = 0") --> Relax{"Relaxation\nchanged dist?"}
    Relax -->|"Yes"| FindCycle{"Find cycle\nin predecessor\ngraph"}
    FindCycle -->|"Cycle found"| Verify{"Is cycle\nnegative?"}
    Verify -->|"Yes"| Yield("Yield cycle 📤")
    Yield --> Relax
    Verify -->|"No"| Relax
    FindCycle -->|"No cycle"| Done("Done ✅")
    Relax -->|"No"| Done
    style Start fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style Relax fill:#fff3e0,stroke:#e65100,color:#2e3440,stroke-width:3px
    style FindCycle fill:#f3e5f5,stroke:#7b1fa2,color:#2e3440,stroke-width:3px
    style Verify fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style Yield fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
    style Done fill:#eceff1,stroke:#455a64,color:#2e3440,stroke-width:3px
</pre>

]

---

### C++ Implementation 💻

Template-based design works with **any graph type**:

```cpp
template <typename DiGraph>
class NegCycleFinder {
    // Graph must support:
    //   for (auto& [node, neighbors] : graph)
    //   for (auto& [target, edge_data] : neighbors)

auto howard(Mapping& dist, Callable get_weight)
        -> cppcoro::generator<Cycle> {
        // Policy iteration loop
        while (relax(dist, get_weight)) {
            for (auto vtx : find_cycle()) {
                co_yield cycle_list(vtx);
            }
        }
    }
};
```

**Graph types**: `std::unordered_map`, `std::list`, custom adapters — any container that yields `(node, neighbors)` pairs.

---

## 📐 Minimum Cycle Ratio (MCR)

---

### Minimum Cycle Ratio Problem 🧮

For a cycle ${\color{lime}C}$ with edge costs ${\color{coral}c}(e)$ and times ${\color{coral}t}(e)$:

$$ \text{ratio}({\color{lime}C}) = \frac{\sum_{e \in {\color{lime}C}} {\color{coral}c}(e)}{\sum_{e \in {\color{lime}C}} {\color{coral}t}(e)} $$

**Goal**: Find the cycle with the **minimum** cost-to-time ratio.

<pre>
graph LR
    A[("a")] -->|"cost=2\ntime=1"| B[("b")]
    B -->|"cost=3\ntime=2"| C[("c")]
    C -->|"cost=1\ntime=1"| A
    style A fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style B fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style C fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
</pre>

]

Ratio of this cycle: $(2 + 3 + 1) / (1 + 2 + 1) = 6 / 4 = 1.5$

---

### Solving MCR via Parametric Search 🔎

Convert to a **parametric feasibility problem**:

For a given ratio ${\color{coral}r}$, define adjusted weights:
$$ {\color{coral}w}_{\color{coral}r}(e) = {\color{coral}c}(e) - {\color{coral}r} \cdot {\color{coral}t}(e) $$

The optimal ratio ${\color{coral}r}^*$ is the smallest ${\color{coral}r}$ such that **no negative cycle** exists with weights ${\color{coral}w}_{\color{coral}r}$.

$$ {\color{coral}r}^* = \min_{{\color{lime}C} \in \text{cycles}} \frac{\sum {\color{coral}c}(e)}{\sum {\color{coral}t}(e)} $$

---

### Algorithmic Visualization 🎯

<pre>
graph LR
    Guess("Guess r") --> Weights("Compute\nwᵣ(e) = c(e) - r·t(e)")
    Weights --> NCF{"Negative\nCycle?"}
    NCF -->|"Yes → r is too low"| Raise("Increase r")
    NCF -->|"No → r is feasible"| Lower("Decrease r")
    Raise --> Weights
    Lower --> Weights
    style Guess fill:#f3e5f5,stroke:#7b1fa2,color:#2e3440,stroke-width:3px
    style Weights fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style NCF fill:#fff3e0,stroke:#e65100,color:#2e3440,stroke-width:3px
    style Raise fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style Lower fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
</pre>

]

---

### C++ MCR Implementation 🚀

```cpp
template <typename Graph, typename T, typename Fn1, typename Fn2, typename Mapping>
auto min_cycle_ratio(
    const Graph& gra, T& r0,
    Fn1&& get_cost, Fn2&& get_time,
    Mapping&& dist, size_t max_iters = 1000)
{
    auto calc_weight = [&](const T& r, const auto& edge) -> T {
        return get_cost(edge) - r * get_time(edge);
    };

auto calc_ratio = [&](const auto& C) -> T {
        return T(total_cost(C)) / total_time(C);
    };

return max_parametric(gra, r0,
        std::move(calc_weight), std::move(calc_ratio),
        std::forward<Mapping>(dist), max_iters);
}
```

Uses `max_parametric` — a general parametric search engine powered by Howard's NCF.

---

## 🔄 Parametric Network Problems

---

### Maximum Parametric Problem 📈

**Formulation**: Maximize ${\color{coral}r}$ such that:

$$ \text{dist}[{\color{firebrick}v}] - \text{dist}[{\color{firebrick}u}] \le \text{distance}(e, {\color{coral}r}) \quad \forall e(u,v) \in {\color{lime}E} $$

<pre>
graph LR
    A[("a")] -->|"d(5,r)"| B[("b")]
    A -->|"d(2,r)"| C[("c")]
    B -->|"d(3,r)"| D[("d")]
    C -->|"d(4,r)"| D
    style A fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style B fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style C fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
    style D fill:#fff9c4,stroke:#f57f17,color:#2e3440,stroke-width:3px
</pre>

]

Where $\text{distance}(e, {\color{coral}r})$ is a **monotone decreasing** function of ${\color{coral}r}$.

---

### Parametric Algorithm 🔧

.font-sm.mb-xs[
```cpp
template <typename Graph, typename T, typename Fn1, typename Fn2, typename Mapping>
auto max_parametric(const Graph& gra, T& r_opt,
    Fn1&& distrance, Fn2&& zero_cancel,
    Mapping&& dist, size_t max_iters = 1000)
{
    auto ncf = NegCycleFinder<Graph>(gra);

for (auto niter = 0U; niter != max_iters; ++niter) {
        // Compute edge weights for current r
        auto get_weight = [&](const auto& edge) {
            return distrance(r_opt, edge);
        };

// Find negative cycle
        for (auto&& c_min : ncf.howard(dist, get_weight)) {
            // Update r_opt from the found cycle
            auto r_min = zero_cancel(c_min);
            if (r_min >= r_opt) return c_opt;
            r_opt = r_min;
            // Update distances along cycle
            for (auto&& edge : c_min) {
                dist[edge.target] = dist[edge.source] - get_weight(edge);
            }
        }
    }
}
```
]

---

## 🧠 Network Oracle Architecture

---

### What is a Separation Oracle? 🎯

In **cutting-plane methods**, an oracle answers:

> "Is point ${\color{green}x}$ feasible? If not, give me a **cutting plane** to eliminate it."

$$ \text{find } {\color{green}x}, {\color{firebrick}u} \quad \text{s.t.} \quad {\color{firebrick}u_j} - {\color{firebrick}u_i} \le h(e, {\color{green}x}), \ \forall e(i,j) \in {\color{lime}E} $$

<pre>
graph LR
    Input("Current point x") --> Oracle("NetworkOracle")
    Oracle --> Decision{"Negative\nCycle?"}
    Decision -->|"No → Feasible ✅"| Done("Return empty")
    Decision -->|"Yes → Infeasible ❌"| Cut("Build cutting plane\n(gradient, intercept)")
    Cut --> Return("Return cut\nfor ellipsoid update")
    style Input fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style Oracle fill:#f3e5f5,stroke:#7b1fa2,color:#2e3440,stroke-width:3px
    style Decision fill:#fff3e0,stroke:#e65100,color:#2e3440,stroke-width:3px
    style Done fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
    style Cut fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style Return fill:#fce4ec,stroke:#ad1457,color:#2e3440,stroke-width:3px
</pre>

]

---

### NetworkOracle C++ Implementation 💪

.font-sm.mb-xs[
```cpp
template <typename Graph, typename Mapping, typename Fn>
class NetworkOracle {
    NegCycleFinder<Graph> _S;  // Howard's algorithm
    Fn _h;                      // Constraint function

template <typename Arr>
    auto assess_feas(const Arr& xval)
        -> std::optional<std::pair<Arr, double>>
    {
        auto get_weight = [&](const auto& edge) {
            return this->_h.eval(edge, xval);
        };

// Find first negative cycle
        for (auto&& C : this->_S.howard(this->_u, get_weight)) {
            auto grad = zeros(xval);
            auto fval = 0.0;
            for (auto&& edge : C) {
                fval -= this->_h.eval(edge, xval);
                grad -= this->_h.grad(edge, xval);
            }
            return {{std::move(grad), fval}};  // Cutting plane!
        }
        return {};  // Feasible
    }
};
```
]

The **cut** is: ${\color{green}g}^T ({\color{green}x}' - {\color{green}x}) + f \le 0$ where ${\color{green}g}$ is the gradient and $f$ the intercept.

---

## ✨ Optimal Matrix Scaling

---

### Problem Statement 📋

Given a matrix ${\color{green}A} = ({\color{coral}a}_{ij})$, find positive diagonal scaling factors ${\color{firebrick}u}_i$ to minimize:

$$ \min \frac{{\color{coral}\pi}}{{\color{coral}\psi}} $$

$$ \text{s.t.} \quad {\color{coral}\psi} \le {\color{firebrick}u}_i \cdot |{\color{coral}a}_{ij}| \cdot {\color{firebrick}u}_j^{-1} \le {\color{coral}\pi}, \ \forall {\color{coral}a}_{ij} \ne 0 $$

<pre>
graph LR
    A[("v₁")] -->|"|a₁₂|"| B[("v₂")]
    B -->|"|a₂₁|"| A
    A -->|"|a₁₃|"| C[("v₃")]
    C -->|"|a₃₁|"| A
    B -->|"|a₂₃|"| C
    C -->|"|a₃₂|"| B
    style A fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style B fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style C fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
</pre>

]

**Orlin & Rothblum (1985)** — reduces to a parametric network problem! 🎯

---

### Log-Scale Transformation 📊

Take logs: ${\color{green}x} = (\log {\color{coral}\pi}, \log {\color{coral}\psi})$, $\text{cost} = \log |{\color{coral}a}_{ij}|$

The constraint becomes:

$$ \log {\color{coral}\psi} \le \log {\color{firebrick}u}_i + \log |{\color{coral}a}_{ij}| - \log {\color{firebrick}u}_j \le \log {\color{coral}\pi} $$

Or equivalently, for each edge $(i,j)$ with cost pair $({\color{coral}a}_{ij}, {\color{coral}a}_{ji})$:

$$ \min({\color{coral}\pi} - {\color{coral}a}_{ji},\ {\color{coral}a}_{ij} - {\color{coral}\psi}) \ge 0 $$

Where $( {\color{coral}\pi}, {\color{coral}\psi} ) = ({\color{green}x_0}, {\color{green}x_1})$ in log scale.

---

### Ratio Oracle ⚖️

```cpp
class Ratio {
    Fn _get_cost;  // extracts (a_ij, a_ji) from edge data

template <typename Edge>
    auto eval(const Edge& edge, const Vec& x) const -> double {
        const auto [aij, aji] = this->_get_cost(edge);
        return std::min(x[0] - aji, aij - x[1]);
    }

template <typename Edge>
    auto grad(const Edge& edge, const Vec& x) const -> Vec {
        const auto [aij, aji] = this->_get_cost(edge);
        if (x[0] - aji < aij - x[1])
            return Vec{1.0, 0.0};   // π-bound active
        return Vec{0.0, -1.0};       // ψ-bound active
    }
};
```

The subgradient tells us **which bound is tighter** — the maximum or minimum ratio.

---

### Complete Pipeline 🏗️

```cpp
// 1. Build graph from matrix sparsity
CostGraph gra = { /* edges with (log|a_ij|, log|a_ji|) */ };

// 2. Create oracle
auto get_cost = [](const auto& edge_data) { return edge_data; };
OptScalingOracle oracle(gra, dist, get_cost);

// 3. Run cutting-plane optimization
Ell ellip(initial_size, x_init);
auto gamma = inf;
auto [xbest, niter] = cutting_plane_optim(oracle, ellip, gamma);
```

<pre>
graph LR
    A("Matrix A") --> B("Graph\n(sparsity)")
    B --> C("OptScaling\nOracle")
    C --> D("Cutting-Plane\n(Ellipsoid)")
    D --> E("Optimal\nScaling ✅")
    style A fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style B fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style C fill:#f3e5f5,stroke:#7b1fa2,color:#2e3440,stroke-width:3px
    style D fill:#fff3e0,stroke:#e65100,color:#2e3440,stroke-width:3px
    style E fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
</pre>

]

---

## 🏭 Applications

---

### Application 1: Clock Skew Scheduling ⏰

In VLSI design, **clock skew** is the difference in clock arrival times between sequential elements.

$$ {\color{firebrick}t}(j) - {\color{firebrick}t}(i) \le T_\text{clk} - {\color{coral}D}(i,j) $$

$$ {\color{firebrick}t}(i) - {\color{firebrick}t}(j) \le {\color{coral}d}(i,j) - T_\text{hold} $$

<pre>
graph LR
    FF1["FF₁\n(clk at t₁)"] -->|"D_max"| FF2["FF₂\n(clk at t₂)"]
    subgraph Constraints["Timing Constraints"]
        C1("t₂ - t₁ ≤ T_clk - D_max")
        C2("t₁ - t₂ ≤ D_min - T_hold")
    end
    FF1 --> Constraints
    FF2 --> Constraints
    style FF1 fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style FF2 fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style Constraints fill:#e3f2fd,stroke:#1565c0,stroke-width:3px
</pre>

]

This is a **difference constraints system** — solvable via negative cycle detection!

---

### Application 2: Performance Analysis 📊

**Minimum cycle ratio** determines the maximum clock frequency of a circuit:

$$ F_{\max} = \frac{1}{\max\text{-cycle-ratio}} $$

Find the cycle where:
$$ {\color{coral}r} = \max_{{\color{lime}C}} \frac{\sum \text{gate delays}}{\sum \text{registers}} $$

This is the **critical cycle** that limits performance. 🏎️

]
.pull-right[
.mermaid[

<pre>
graph LR
    A["R₁"] -->|"2ns"| B["R₂"]
    B -->|"3ns"| C["R₃"]
    C -->|"1ns"| D["R₄"]
    D -->|"2ns"| A
    B -->|"1ns"| D
    style A fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style B fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style C fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
    style D fill:#fff9c4,stroke:#f57f17,color:#2e3440,stroke-width:3px
</pre>

]

Cycle R₁→R₂→R₃→R₄→R₁ has ratio $(2+3+1+2)/(4) = 2.0\text{ns}$
]

---

### Application 3: Matrix Balancing ⚖️

**Optimal matrix scaling** has applications in:

- **Preconditioning**: Improve numerical stability of matrix computations
- **Economic equilibrium**: Find exchange rates that balance trade
- **Traffic networks**: Scale origin-destination matrices
- **Chemical engineering**: Balance reaction rates

<pre>
graph LR
    subgraph Before["Before Scaling"]
        M1["|a₁₁|  |a₁₂|  |a₁₃|\n|a₂₁|  |a₂₂|  |a₂₃|\n|a₃₁|  |a₃₂|  |a₃₃|"]
    end
    subgraph After["After Scaling"]
        M2["u₁·|a₁₁|·u₁⁻¹  u₁·|a₁₂|·u₂⁻¹  u₁·|a₁₃|·u₃⁻¹\nu₂·|a₂₁|·u₁⁻¹  u₂·|a₂₂|·u₂⁻¹  u₂·|a₂₃|·u₃⁻¹\nu₃·|a₃₁|·u₁⁻¹  u₃·|a₃₂|·u₂⁻¹  u₃·|a₃₃|·u₃⁻¹"]
    end
    Before -->|"Find uᵢ\nto minimize\nπ/ψ"| After
    style Before fill:#fff3e0,stroke:#e65100,stroke-width:3px
    style After fill:#e8f5e9,stroke:#2e7d32,stroke-width:3px
</pre>

]

---

## 📦 Project Structure

---

### Repository Layout 📁

```
digraphx-cpp/
├── include/digraphx/
│   ├── neg_cycle.hpp        # Howard's algorithm
│   ├── parametric.hpp       # Parametric search
│   ├── min_cycle_ratio.hpp  # MCR algorithm
│   ├── min_parametric_q.hpp # Discrete parametric
│   └── neg_cycle_q.hpp      # Queue-based NCF

netoptim-cpp/
├── include/netoptim/
│   ├── network_oracle.hpp    # Separation oracle
│   ├── optscaling_oracle.hpp # Matrix scaling oracle
│   ├── neg_cycle.hpp         # Local NCF (legacy)
│   ├── parametric.hpp        # Parametric (legacy)
│   ├── min_cycle_ratio.hpp   # MCR (legacy)
│   └── primal_dual.hpp       # Vertex cover
└── test/source/
    └── test_*.cpp            # Test suite
```

---

### Build & Test 🔨

**xmake** build system:

```bash
# Build and run all tests
xmake -y -j 10
xmake run test_netoptim

# Run specific test
xmake run test_netoptim -tc="Test *"
```

**17 test cases** covering:
- ✅ Negative cycle detection
- ✅ Minimum cycle ratio
- ✅ NetworkOracle separation
- ✅ Optimal matrix scaling
- ✅ Large graph stress tests

All passing with **Howard's algorithm** from `digraphx-cpp`! 🎉

---

### Dependency Graph 🔗

<pre>
graph TD
    NETOPTIM["netoptim-cpp\n🏗️ Network Oracles"] --> DIGRAPHX["digraphx-cpp\n🔄 NegCycleFinder"]
    NETOPTIM --> ELLALGO["ellalgo-cpp\n🎯 Ellipsoid Method"]
    DIGRAPHX --> PY2CPP["py2cpp\n📦 Python-like containers"]
    DIGRAPHX --> CPPCORO["cppcoro\n🔄 Coroutine generators"]
    NETOPTIM --> XNETWORK["xnetwork-cpp\n🌐 Graph structures"]
    style NETOPTIM fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style DIGRAPHX fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style ELLALGO fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
    style PY2CPP fill:#fff9c4,stroke:#f57f17,color:#2e3440,stroke-width:3px
    style CPPCORO fill:#f3e5f5,stroke:#7b1fa2,color:#2e3440,stroke-width:3px
    style XNETWORK fill:#e0f2f1,stroke:#00695c,color:#2e3440,stroke-width:3px
</pre>

]

---

## 📈 Performance & Benchmarks

---

### Why Howard's Algorithm? 🏅

| Algorithm | Source Needed | Memory | Best For |
|-----------|:---:|:------:|----------|
| **Bellman-Ford** | ✅ Yes | $O(V)$ | Sparse graphs |
| **SPFA** | ✅ Yes | $O(V)$ | Average-case speed |
| **Howard (Policy Iteration)** | ❌ No | $O(V)$ | **Dense graphs, many cycles** |

**Key advantage**: Howard's algorithm **reuses distance estimates** between iterations — no need to restart from scratch. The generator pattern yields cycles on demand without storing all of them. 🚀

---

### Complexity Analysis 📊

- **Howard's NCF**: $O(V \cdot E \cdot C)$ where $C$ = number of policy iterations
  - Typically converges in very few iterations
- **Parametric search**: $O(\text{iters} \times \text{NCF complexity})$
  - Usually 5–20 iterations in practice
- **Cutting-plane method**: $O(n^2 \times \text{oracle cost})$ for ellipsoid method

<pre>
graph TD
    subgraph Cost["Cost Breakdown"]
        A("Parametric Search\n5-20 iterations") --> B("NCF per iteration\nO(V·E·C)")
        B --> C("Total O(V·E·C·iters)")
    end
    style Cost fill:#e3f2fd,stroke:#1565c0,stroke-width:3px
    style A fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style B fill:#fff3e0,stroke:#e65100,color:#2e3440,stroke-width:3px
    style C fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
</pre>

]

---

## 🔮 Future Directions

---

### What's Next? 🚀

#### Algorithmic
- **Negative Cycle Finder Q** — queue-based variant for sparse graphs
- **Discrete parametric** optimization
- **Multi-parameter** network problems
- **Planar graph** specializations

]

#### Implementation
- **Python bindings** via pybind11
- **GPU acceleration** for large matrices
- **WebAssembly** port for browser demos
- **Better benchmarks** & visualization

]

---

### Key Takeaways 🎯

1. **Negative cycles** are the fundamental primitive for network optimization
2. **Howard's algorithm** provides efficient, generator-based detection
3. **Parametric search** reduces ratio problems to feasibility checks
4. **Cutting-plane oracles** bridge network optimization and convex optimization
5. **Optimal matrix scaling** is a beautiful application of this framework

> "All models are wrong, but some are useful." — George Box

The network optimization framework in `digraphx-cpp` and `netoptim-cpp` provides **powerful, reusable primitives** for solving a wide range of practical problems! 🌊

---

## 🚀 New: Fast Variants

---

### digraphx-fast & netoptim-fast ⚡

Two new projects optimizing for **large-scale** network optimization:

.font-sm.mb-xs[

| | Original | Fast |
|---|---|---|
| **Repository** | `digraphx-cpp` / `netoptim-cpp` | `digraphx-fast` / `netoptim-fast` |
| **Graph storage** | Hash-based (unordered_map) | **CSR format** (contiguous arrays) |
| **Cycle enumeration** | `cppcoro::generator` (coroutines) | Callback-based (`Fn&&`) |
| **Pred tracking** | `unordered_map<Node, pair<Node, Edge>>` | `vector<node_t>` + `vector<size_t>` |
| **Visited tracking** | `unordered_map<Node, Node>` (heap alloc) | `vector<node_t>` (start-vertex ID, no alloc) |
| **Warm-start** | ❌ No | ✅ `howard_warm()` preserves predecessor |
| **max_cycles** | ❌ Yields all cycles | ✅ Configurable (default 1) |
| **Dependencies** | cppcoro, py2cpp, xnetwork | **None** (STL only) |

]

---

### CSR Graph Storage 📦

**Before** (hash-based): each node has a `unordered_map` — per-node allocation, cache misses

.font-sm.mb-xs[
```cpp
// Old: unordered_map per node, pointer-chasing iteration
std::unordered_map<uint32_t, std::unordered_map<uint32_t, double>> graph;
```
]

**After** (CSR): single flat arrays — linear scan, cache-friendly

.font-sm.mb-xs[
```cpp
// New: three contiguous vectors
struct CSRGraph {
    std::vector<uint32_t> offsets;   // [V+1]
    std::vector<uint32_t> targets;   // [E]
    std::vector<double>   weights;   // [E]
};
```
]

Relaxation becomes a tight linear loop:

.font-sm.mb-xs[
```cpp
for (node_t u = 0; u < V; ++u) {
    auto du = dist[u];
    for (size_t e = offsets[u]; e < offsets[u + 1]; ++e) {
        auto d = du + weights[e];
        if (dist[targets[e]] > d) { /* update */ }
    }
}
```
]

---

### Epoch-Free Visited Tracking 👁️

The old `_find_cycle` used `std::unordered_map<Node, Node>` — **heap allocated every call**.

The fast version uses a **pre-allocated array** storing the originating start-vertex:

```cpp
// O(V) memset instead of O(V) unordered_map inserts
auto _find_cycle() -> std::optional<node_t> {
    std::fill(_visited.begin(), _visited.end(), node_t(-1));
    for (node_t v = 0; v < V; ++v) {
        if (_visited[v] != node_t(-1)) continue;
        auto u = v;
        _visited[u] = v;  // mark: "visited by start v"
        while (_pred_node[u] != node_t(-1)) {
            u = _pred_node[u];
            if (_visited[u] == v) return u;   // cycle back to v
            if (_visited[u] != node_t(-1)) break;  // visited by other start
            _visited[u] = v;
        }
    }
    return {};
}
```

No hash table, no allocations, no `clear()` — just a single `std::fill`.

---

### Multi-Edge Support 🔄

**Both versions fully support parallel edges** (multiple edges between the same node pair).

.pull-left[
In the old version, the neighbor container (e.g. `std::list`) naturally holds duplicates:

.font-sm.mb-xs[

```cpp
// Old: two parallel edges 2→0
{2, {{1, 1.0}, {0, 2.0}, {0, 1.0}}}
//                     ^^^^^^^^  ^^^^^^^^  parallel!
```
]
]
.pull-right[
In the new CSR version, each `add_edge()` creates a separate CSR entry:

.font-sm.mb-xs[
```cpp
builder.add_edge(2, 0, 2.0);  // edge index 3
builder.add_edge(2, 0, 1.0);  // edge index 4

```
]

]

**Advantage of the new version**: cycles return **edge indices** (not edge data values), so the caller can disambiguate which specific parallel edge participates in the cycle.

<pre>
graph LR
    A[("2")] -->|"e₃: 2.0"| B[("0")]
    A -->|"e₄: 1.0"| B
    style A fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style B fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
</pre>

]

During relaxation, both edges are evaluated. The lower-weight edge wins (`_pred_edge[0] = 4`), correctly selecting the shortest path.

---

### Performance Trade-Offs ⚖️

| Criterion | Pick Original | Pick Fast |
|---|---|---|
| Graph type flexibility | ✅ Any iterable | ❌ CSR only |
| Coroutine generators | ✅ `cppcoro::generator` | ❌ Callback API |
| Multi-edge disambiguation | ⚠️ By value | ✅ By edge index |
| Large graphs (1M+ edges) | ❌ Hash overhead | ✅ CSR cache-friendly |
| Memory constrained | ❌ Per-node hash tables | ✅ Single allocation |
| Parametric search warm-start | ❌ Restarts each iteration | ✅ Reuses predecessor |
| Streaming / dynamic graphs | ✅ Easy modification | ⚠️ Needs rebuild |

---

### New Dependency Graph 🔗

<pre>
graph TD
    subgraph Fast["⚡ Fast Variants"]
        FAST_NET["netoptim-fast\n🏗️ Network Oracles"] --> FAST_DIGRAPH["digraphx-fast\n🔄 CSR NegCycleFinder"]
        FAST_NET -.->|"runtime"| ELLALGO["ellalgo-cpp\n🎯 Ellipsoid Method"]
    end
    subgraph Original["📦 Original"]
        NETOPTIM["netoptim-cpp\n🏗️ Network Oracles"] --> DIGRAPHX["digraphx-cpp\n🔄 NegCycleFinder (generator)"]
        DIGRAPHX --> CPPCORO["cppcoro\n🔄 Coroutines"]
    end
    FAST_DIGRAPH -->|"no deps"| STL["C++20 STL only"]
    NETOPTIM --> ELLALGO
    DIGRAPHX --> PY2CPP["py2cpp\n📦 Containers"]
    NETOPTIM --> XNETWORK["xnetwork-cpp\n🌐 Graph structures"]
    style FAST_NET fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style FAST_DIGRAPH fill:#ffcdd2,stroke:#c62828,color:#2e3440,stroke-width:3px
    style NETOPTIM fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style DIGRAPHX fill:#e3f2fd,stroke:#1565c0,color:#2e3440,stroke-width:3px
    style ELLALGO fill:#c8e6c9,stroke:#2e7d32,color:#2e3440,stroke-width:3px
    style STL fill:#f3e5f5,stroke:#7b1fa2,color:#2e3440,stroke-width:3px
</pre>

]

> **Note**: `digraphx-fast` is truly STL-only (no external dependencies). `netoptim-fast` headers are also STL-only — they implement the oracle interface. But to *solve* optimization problems you plug them into `ellalgo-cpp`'s `cutting_plane_optim()` / `Ell` at runtime (dashed edge).

---

# 🎤 Q&A

## Thank You! 🙏

@luk036 👨‍💻 · 2026 📅

**References**:
- Orlin & Rothblum (1985) — Optimal Matrix Scaling
- Howard (1960) — Dynamic Programming and Markov Processes
- Karp (1967) — Minimum Mean Cycle Algorithm