I implemented contacts warm starting the way Erin Catto describes in his 2009 GDC Box2D lecture (box2d.org/files/GDC2009/GDC2009_Catto_Erin_Solver.ppt).

However this seems to only work for inelastic contacts (no bouncing).

Can anyone point me in the right direction for how to extend this idea to incorporate elastic collisions?

I tried looking through the Box2D source for this, but I could not figure out how it works. Any help is appreciated!

Ok. I found the answer in the source code for Irlan Robson's Bounce Physics engine.

https://github.com/irlanrobson/bounce/blob/master/src/bounce/dynamics/contacts/contact_solver.cpp

The key is to store the restitution velocity bias before the solver starts (for example in the contacts setup method for the solver).

if (relativeNormalVelocity < -THRESHOLD) {
	contact.velocityBias = -relativeNormalVelocity * restitutionCoefficient;
} else {
	contact.velocityBias = 0;
}

Then in each step of the solver when you calculate the accumulated impulse you add in the velocity bias:

float oldImpulse = contact.impulse;

float impulse = (-relativeNormalVelocity + contact.velocityBias) * contact.mass;
contact.impulse = max(contact.impulse + impulse, 0.0f);

impulse = contact.impulse - oldImpulse;

That way the relative velocity for the restitution is preserved throughout the iterations.

Thank you @Irlan Robson