Brian Jackson

This is a good question, hopefully my explanations below provide some insight: ```julia ForwardDiff.jacobian(testfoo1,q1) # equivalent to ForwardDiff.jacobian(normalize, q1) ForwardDiff.jacobian(testfoo1,q2) # 9-dimensional input since all 3D Rotations are `StaticMatrix{3,3}`, so...

At least in my fields of dynamics, aerospace, robotics, and mechanical engineering, it's more standard to have the quaternion ordered as `[w,x,y,z]`. So I'm in favor of keeping it as-is....

It happens with both input and output archives. I can certainly use a `std::span` for output, but what about for input? Is there a way to customize the behavior to...

 Found that most of the time is spent in forming Ahat and Bhat. Spawned #23.

 New results after #26. Next steps are to further optimize `updateBhat!` by maybe doing some more caching, and further optimizing `eval_c!`. It seems a lot of time is being...

It looks like the penalty parameter needs to be sufficiently high in order for it to converge.

Not clear how exactly to implement this in the bilinear case. May look into the motivation behind it to see if it can be adapted to the bilinear case.