A Computer Scientist Walks Into a Math Department

Divide and conquer, applied to an interval instead of an array

In 2005 I had six months and nothing to do with them.

My admission to Cornell had been deferred to the spring, I was 25, and I was living in Lahore with a gap in my life the exact shape of a semester. The sensible move would have been to take a contract job and bank the money. Instead I walked across campus at LUMS, into the mathematics department, and asked a professor whether he needed anyone to write code.

That decision produced the only academic paper I've ever published — and, more interesting to me now, it produced a lesson about what happens when you drop a computer scientist into a problem that mathematicians have been staring at.

The code I wrote for it survived. I found the archive recently, twenty years on, and I've put it on GitHub along with a modern C++20 port. Reading it again has been humbling in about four different ways, and I'll get to all of them.

TL;DR

• A math professor had a genuinely elegant method for solving stiff ODEs — Sobolev gradients — that avoids time-stepping entirely.
• It worked beautifully on short intervals and fell apart on long ones. On u' = u, going from [0,1] to [0,5] took the solver from 57 steps to 148,798.
• While implementing it, I noticed the failure had the shape of a problem I already knew: divide and conquer. Split the interval, recurse, and — the part I'm still proud of — use the cost of the subproblem itself as the splitting criterion.
• That turned a method that needed hand-tuning into one that derives its own structure from the problem. [0,5] came down from 148,798 steps to 487.
• The paper was published in the International Journal of Computer Mathematics in 2008. The mathematics wasn't mine; the contribution I brought to it came from CS.
• Porting the code in 2026 uncovered a sign error in my conjugate gradient solver that has been sitting inside the published tables for eighteen years, quietly explaining a number nobody questioned.

The problem: equations that punish you for marching

If you want to solve a differential equation numerically, the obvious approach is to march. Start at a known point, take a small step forward, compute the next value, repeat. Euler's method is this in its simplest form, and Runge–Kutta is this done carefully.

Marching works right up until the equation is stiff.

Stiffness is one of those concepts with no single clean definition, but the practical version is: the equation contains two or more wildly different time scales. Something in the system varies violently fast, while the solution you actually care about is nearly flat.

The canonical example — the one Cleve Moler uses to explain stiffness in MATLAB, and the one we ended up testing against — is the flame-size equation:

du/dt = u² - u³,    u(0) = 10⁻⁴

This models a match head igniting. The flame grows slowly, slowly, slowly, and then at around t = 10⁴ it jumps to a stable size and flattens out forever. Here's the cruelty: an explicit solver has to take steps small enough to survive the jump the entire time, including across the ten thousand units of nearly nothing that precede it. The step size is dictated by the worst moment in the system, not the current one.

You can pay for that with implicit methods, which is what MATLAB's ode15s does. You solve an equation at every step instead of just evaluating one. It works. It's also the standard answer, and standard answers are less fun.

The professor's idea: stop marching entirely

Sultan Sial's approach — building on Sobolev gradient theory developed by J.W. Neuberger, who literally wrote the book on it — was to refuse the premise.

Don't march. Solve the whole interval at once.

Write the ODE as L(u) = 0, where u is now not a value at a point but a vector holding the entire solution across the interval. Then define an error functional:

F(u) = ½ ⟨L(u), L(u)⟩

F(u) is zero exactly when u solves the equation everywhere. So finding the solution becomes minimizing a function — and we have a well-worn tool for that. Steepest descent. Start with a guess, compute the gradient, step downhill, repeat.

And here's the part I find genuinely beautiful: there is no time-stepping anywhere in this. You're not evolving forward through time. You're taking a wrong answer covering the whole interval and pushing it, everywhere at once, toward being right. No time-stepping means no notion of a time scale, and stiffness is a statement about time scales. The thing that makes the problem hard has nowhere to attack.

The catch is which downhill you follow

"Compute the gradient" hides the entire difficulty, because a gradient is not an absolute thing. A gradient is the direction of steepest increase per unit change in the argument — and "unit change" depends on how you measure size. Change the inner product and you change what "steepest" means.

Do the obvious thing — descend in L² — and you walk straight into the CFL condition. Refine your grid and the step size has to collapse. You've smuggled the same instability back in through the side door.

Instead, measure with an inner product that accounts for the derivative too:

⟨u,v⟩_S = ⟨D₀u, D₀v⟩ + ⟨D₁u, D₁v⟩

where D₀ averages adjacent nodes and D₁ differences them. The gradient with respect to that metric — the Sobolev gradient — is dramatically smoother, and the step size stops caring how fine your grid is.

The price is that you don't get it for free. Getting the Sobolev gradient from the L² one means solving a linear system at every single descent step:

(D₀ᵀD₀ + D₁ᵀD₁) ∇_S F = ∇_{L²} F

which is what the conjugate gradient solver in my code is for. Hold that thought — it comes back to bite me in about fifteen paragraphs.

The payoff is not subtle. Same problem, same grid, two different notions of "downhill":

Published results, Table 1. Note the direction of the Sobolev column as the grid refines.

Look at the last column again. As the grid gets four times finer, the L² method needs 13× more steps, and the Sobolev method needs fewer. That's the CFL condition being absent rather than merely mitigated.

Where it fell apart

So the method is elegant and the numbers are excellent. Where's the problem?

Lengthen the interval.

Keeping the node spacing fixed at 0.01 and just solving u' = u over longer and longer stretches:

Published results, Table 2. This is not a gentle degradation.

On a log scale the bars step up almost evenly — that's what constant multiplicative growth looks like

Each unit of interval costs roughly 6× the last. That is a method with a hard practical ceiling, and "stiff problems on long intervals" is precisely the use case you'd want it for. The flame-size equation lives on [0, 20000].

Three things go wrong at once, and only the third is interesting.

The first two are boring: the linear solve gets more expensive on bigger arrays, and the forward-Euler starting guess gets worse the further it has to run. Real, but unsurprising.

The third one is the actual disease. The gradient is local.

Picture the solver halfway through. The left half of your interval has essentially converged — those values are right. The right half is still garbage. Now compute the gradient. It's near zero across the converged left half, near zero across the not-yet-organized right half, and it spikes at the seam where they meet.

So descent does almost nothing useful. Worse: the correct left half gets dragged around by the disturbance at the seam while the right half slowly sorts itself out. You are actively un-solving the part you already solved. The bigger the interval, the more of the run is spent in this state.

The L² gradient is zero except at the seam. The Sobolev gradient — the smooth one that makes the whole method work — rises across the entire half that was already correct

There's an irony in that bottom trace worth sitting with. Smoothing the gradient is exactly what defeats the CFL condition and makes this method viable. It's also precisely why the converged half gets disturbed: a smooth gradient spreads the correction, and spreading it means spreading it onto places that didn't need correcting. The property that saves the method on short intervals is the one that sinks it on long ones.

The part where I was useful

Here is where I want to be honest about what my contribution actually was, because it wasn't mathematics.

I was implementing this. I'd been staring at the degradation table for days, and what I kept thinking was: this looks like a problem I know. Not from numerical analysis — from algorithms class. A big instance is expensive and behaves badly. Small instances are cheap and behave well. The expensive thing is made of smaller things.

That's mergesort. That's quicksort. That's binary search. That's the first three weeks of any algorithms course.

The professors already knew that solving on subintervals helped — solve [0,1], take the endpoint value, use it as the initial condition for [1,2], and so on. Doing that with fixed subintervals of length 1 brought [0,5] from 148,798 steps down to 364. Problem solved, apparently.

Except it isn't, and this is the bit I care about. It just moves the question. How long should the subintervals be?

Length 1 works for u' = u. It's meaningless for the flame equation, where the interval is 20,000 long and nothing happens for the first half of it. And even within a single problem the right answer varies from place to place — a flat stretch wants enormous subintervals, a sharp transition wants tiny ones. You're back to hand-tuning a parameter per problem, by feel, which is exactly the kind of thing that makes a method impressive in a paper and useless in practice.

So the contribution is this: stop choosing.

Don't pick a subinterval length. Try to solve the whole thing. If it costs more than a fixed budget of work without converging, cut it in half and recurse — left first, feeding the left solution's endpoint in as the right half's initial condition.

function sobolevRecursive(interval, maxEffort) returns solution
    solution = solve(interval, maxEffort)
    if (solution is complete)
        return solution
    if (solution is incomplete) or (solution is spiralling out)
        split interval => leftInterval, rightInterval
        leftSolution  = sobolevRecursive(leftInterval,  maxEffort)
        rightSolution = sobolevRecursive(rightInterval, maxEffort)
        solution = combine(leftSolution, rightSolution)
    return solution
end function

Twelve lines. Any CS undergraduate would recognize the shape immediately.

The one design decision worth arguing about is the splitting criterion. Usually you divide and conquer on size — half the array, half the range. Here we divide on cost: work is measured as nodes × minimization steps, and when that exceeds 10,000 without convergence, split. The subproblem's own difficulty decides whether it gets subdivided.

That's what makes it adaptive rather than merely recursive. The recursion goes deep where the problem is hard and stays shallow where it's easy, and nobody has to know in advance where the hard parts are.

Same trick applies to the grid. The formulation was extended to non-uniform node spacings, so nodes get inserted where consecutive values differ too much and removed where the solution is flat.

Published results, Tables 2 and 3.

Read that table carefully, because it's easy to misread as a defeat. Recursive splitting is slightly worse than hand-picked subintervals — 487 against 364. That's the price of not being told the answer. What you buy for that 34% is that nobody had to know length 1 was right.

And on a problem where nobody could have guessed, that trade pays enormously. Turned loose on the flame equation, the algorithm produced this profile on its own:

  start   9,375.00  length 312.50
  start   9,843.50  length  78.00
  start   9,960.50  length   9.50
  start  10,000.00  length   2.25     <- the flame jump
  start  10,009.75  length   4.88
  start  10,078.00  length  78.00
  start  10,312.50  length 312.50

Subintervals 312 units wide across the flat stretch, narrowing to 2.25 across the jump, widening back out afterward. Nobody told it the jump was at 10,000. It found the interesting part of the problem by noticing where the work got expensive.

That's the whole idea in one block of numbers, and I still think it's lovely.

The departments that didn't talk

Here's the institutional part, and I suspect it's the most transferable thing in this post.

At LUMS in 2005, the mathematics department and the computer science department barely spoke. Not out of hostility — just the ordinary drift of two groups with different buildings, different seminars, different conferences, different vocabulary for the same ideas. I was a CS student. Working with a math professor was mildly eccentric.

And the paper that came out of it is entirely a product of both. The mathematics is Sial's and Neuberger's — the functional, the inner product, the whole apparatus for descending in a Sobolev space. None of that is mine, and I couldn't have derived any of it. The recursion is CS, and it came from an undergraduate algorithms reflex.

There was an earlier version of this work, covering the flame-size problem with the standard method and no recursion. What the published paper adds on top of it is the divide-and-conquer scheme — twelve lines of undergraduate computer science, sitting on top of a great deal of mathematics that wasn't mine.

What strikes me now is how close those two ideas were the whole time and how unlikely the collision was. Divide and conquer isn't advanced. It isn't a clever trick I invented. It's first-year material that happened not to be in the room.

I've seen the same pattern in every company I've worked at since. Not math and CS — data science and platform, security and product, the people who know why the query is slow and the people who know why the query exists. The gap is rarely knowledge. It's that nobody's sitting in both rooms.

If you're ever offered six months of nothing, go sit in the other room.

Writing C++ in 2005, without a safety net

Some context on what building this actually looked like, because the working conditions have changed more than the algorithm has.

Stack Overflow did not exist. It launched in 2008, three years later. There was no Copilot, no ChatGPT, no Cursor, no autocomplete worth the name. There was MSDN, there were books, there was Google — and there was reading the code you'd already written until you understood why it was wrong.

I chose C++ for one reason: the professor wanted to compare our implementation's performance against MATLAB's. That's a drag race, and I wanted the biggest engine on the grid. C++ was the V8.

The whole thing is Visual C++ 6, a console project, last modified October 2005. About 2,400 lines. And it's hand-rolled in a way that I don't think I've had to be since:

class DoubleArray
{
public:
    DoubleArray();
    DoubleArray(int size);
    virtual ~DoubleArray() { if(n>0) delete [] ptr; }
    int getSize() { return n; }
    double &operator[](int index) { return ptr[index]; }
    ...

That's a dynamic array. With node insertion and deletion, because the adaptive grid needs it. I wrote my own, because std::vector support in VC6 was not something I trusted for numerical work, and because that's simply what you did.

There is also this, sitting near the top of the header, which made me laugh out loud twenty years later:

#define abs(x) (x<=0?-1.0*x:x)

A macro shadowing a standard function, evaluating its argument twice. Every code review I've run in the last fifteen years would reject that on sight. It works fine here, right up until someone writes abs(i++).

Changing a parameter meant editing main() and recompiling. Selecting which problem to run meant commenting out function calls. There are four near-duplicate variants of the interval-wise solver — sobolevIntervalWise, sobolevIntervalWiseVS, and two overloads — that differ mainly in which of eleven parallel arrays they thread through their signatures. Eleven. Threaded by hand through every call site.

I don't show you this to be self-deprecating. I show you it because it was fast, it was correct enough to publish, and a 25-year-old wrote it alone with no tooling. The constraints were real and the work got done anyway.

Twenty years later, the port

I took the archive and built a modern version: C++20, CMake, std::vector, a real command-line interface, and — for the first time in this code's life — tests.

The port reproduces the paper. 55 minimization steps against a published 57 on the prototypical problem. 3,661 against 3,432 on the flame problem, with matching node counts. The strongest check isn't a table, though — it's that the 2005 binary left a run log sitting in its Debug/ folder:

REC: Spliting at x = 10000
REC: Spliting at x = 4999.5
REC: Spliting at x = 7498.75
REC: Spliting at x = 8747.87
REC: Spliting at x = 9372.94
REC: Spliting at x = 9685.97
REC: Spliting at x = 9841.98

The actual output of the 2005 executable, recovered from the build directory

That's the recursion binary-searching its way toward the flame jump at 10,000. The modern port splits at 9375, 9687.5, 9843.5 — the same cascade, from code rewritten from scratch two decades later. That's the kind of agreement you can't fake.

Then the tests started finding things.

The bug in the published tables

The one that matters is four characters long:

beta = (-1.0*d1)/d0;
d = r + (d*beta);

Textbook conjugate gradient uses beta = +d1/d0. With the sign flipped, the search directions are no longer conjugate — and conjugacy is the entire point of CG. It's what guarantees termination in at most n steps on an n-dimensional system. Flip that sign and you still have a descent method, which is why everything converged and why nobody noticed. You just no longer have conjugate gradient.

Measured on a 40-node system: the corrected version terminates in 42 iterations with a residual of 5e-14. The 2005 version needs more than 100,000 to reach six digits.

And here's the thing that genuinely got me. Go back to the published paper and look at the CG counts: 33,490 CG iterations for 57 minimization steps on a 101-node system. That's roughly 590 CG iterations per solve on a system where CG should need at most 101.

That number was printed in a peer-reviewed journal in 2008. It's anomalous on its face, and it sat there for eighteen years — through my own eyes, two professors, and a review process — because it was in a column nobody was looking at. The paper's argument was about minimization steps. The CG column was housekeeping.

Running it today, on the same problem:

u' = u on [0,5], 501 nodes, recursive mode. Live from the port.

The same answer, for a ninth of the work

Eleven times fewer CG iterations. Six times faster. Same answer.

The conclusions of the paper survive completely — the recursive method still fixes the long-interval problem, the numbers all still point the same direction. But the implementation was leaving a 6× speedup on the table, and the evidence was published.

The port keeps the 2005 behavior as the default, because that's what actually produced the paper. --numerics corrected turns on the fixes. It felt wrong to quietly repair history.

There were others: the recursive function has a code path with no return statement at all, the right subinterval never actually received the left's endpoint the way the paper describes, and the adjoint of D₁ is only correct on a uniform grid — which is fine for the fixed-grid tables and quietly wrong for the non-uniform ones that were the entire point of half the paper. They're all written up in the porting notes.

Four honest defects in code that produced a published result. That's a useful thing to sit with.

The more things change

One last thing, and it's the reason I find this whole exercise more than nostalgia.

I picked C++ in 2005 because I needed to beat MATLAB and nothing else was fast enough. Java existed. Python existed. Neither was going to win a numerical drag race.

Twenty years later, the world runs on Python — and NumPy is C, SciPy is C and Fortran, PyTorch is C++, Polars is Rust. Every time you write a line of Python that has to be fast, you are calling into something that somebody wrote in a systems language, for exactly the reasons I picked C++ in a lab in Lahore.

The languages we type in changed completely. The layer underneath didn't move.

Conclusion

What I keep returning to isn't the algorithm. It's that the entire contribution was noticing a shape.

The mathematics was beyond me — I couldn't have derived the Sobolev gradient formulation, then or now. What I had was a different set of reflexes: an instinct that says expensive problem, cheap subproblems, recurse, and a second instinct that says don't make the user pick the parameter, measure it instead. Those aren't advanced ideas. They're the plumbing of computer science, and they were worth something precisely because they weren't in that room.

Twenty years later, my tests found four bugs in the code that made it work. The idea was right and the implementation was flawed, and both of those things stay true at the same time. I'm oddly comforted by that. It's the most accurate description of engineering I know.

Also, I'd like the record to show that in 2005 I hand-rolled a dynamic array with node insertion, and it still runs.

The code: github.com/danishm/steepest-descent-2005 — the 2005 sources untouched, plus the modern C++20 port with tests and porting notes. MIT licensed.

The paper: D. Mujeeb, J.W. Neuberger and S. Sial, Recursive form of Sobolev gradient method for ODEs on long intervals, International Journal of Computer Mathematics 85:11 (2008), 1727–1740. doi:10.1080/00207160701558465