Fix high-order Trotter bug

[guanghaolow]

Some tests by Jan Gukelberger revealed that the error of the 4th order Trotter formulas are not scaling as expected. The error can be traced to a difference in notation. In the cited reference, the order of Trotter formulas are indexed by 2k. In this library, order is indexed by k.

[jgukelberger]

Thanks for the quick fix!

[cgranade]

Before merging this in, we should make sure we have a documented regression check, ideally in the form of a unit test.

[cgranade]

If you'd like, @guanghaolow, I'm happy to help get regression tests added. Thanks!

[jgukelberger]

@cgranade @guanghaolow - any updates on this? Just coming back to this project, I was a little baffled to find this fix hadn't made it into any release yet. While I appreciate the goal of improving test coverage, it's maybe not ideal if the turnaround time of one-line fixes is blown up from days to months?

Regarding a regression test, I don't have great ideas currently. I found the bug by plotting Trotter error vs step size, where the error is computed with the help of some state dumper hacking and scipy.linalg. At some point that project may become a sample, which would then also be sensitive to the bug, but it's not something that could be quickly turned into a good test at the moment.

[cgranade]

@jgukelberger: I definitely understand where you're coming from, but since this operation is pretty widely used, I want to make sure the fix is correct before we merge it in. In particular, looking back at quant-ph/0508139, it looks as though the definition of 𝑆₂(λ) matches _Trotter2ImplCA:

    operation _Trotter2ImplCA<'T> ((nSteps : Int, op : ((Int, Double, 'T) => Unit is Adj + Ctl)), stepSize : Double, target : 'T) : Unit is Adj + Ctl
    {
        for (idx in 0 .. nSteps - 1)
        {
            op(idx, stepSize * 0.5, target);
        }
        
        for (idx in nSteps - 1 .. -1 .. 0)
        {
            op(idx, stepSize * 0.5, target);
        }
    }

I could be misunderstanding, but that suggests to me that the indexing conventions are the same in the paper and in the library. Similarly, the version before this PR seemed to match the indexing convention for 𝑆₂ₖ(λ) and 𝑝ₖ as well:

let stepSizeOuter = _TrotterStepSize(order);
let stepSizeInner = 1.0 - 4.0 * stepSizeOuter;
_TrotterArbitraryImplCA(order -2, (nSteps, op), stepSizeOuter * stepSize, target);
_TrotterArbitraryImplCA(order -2, (nSteps, op), stepSizeOuter * stepSize, target);
_TrotterArbitraryImplCA(order -2, (nSteps, op), stepSizeInner * stepSize, target);
_TrotterArbitraryImplCA(order -2, (nSteps, op), stepSizeOuter * stepSize, target);
_TrotterArbitraryImplCA(order -2, (nSteps, op), stepSizeOuter * stepSize, target);

function _TrotterStepSize (order: Int) : Double {
    return 1.0 / (4.0 - PowD(4.0, 1.0 / (IntAsDouble(order) - 1.0)));
}

Given that, I was hoping to get @guanghaolow's input to make sure that this is the correct fix, rather than risking a second break. That said, as I commented above I'm also quite happy to do what I can to help get this bug fix in. I just want to be cautious about the fix. Thanks for reporting, and for your understanding!

[guanghaolow]

Thanks for pointing out the difference between the paper quant-ph/0508139 and the proposed change Chris!

We can better match the notation of the paper by putting the missing factor of two into
let stepSizeOuter = _TrotterStepSize(order/2);
&
return 1.0 / (4.0 - PowD(4.0, 1.0 / (IntAsDouble(2*order) - 1.0)));

instead of the proposed change
let stepSizeOuter = _TrotterStepSize(order);
&
return 1.0 / (4.0 - PowD(4.0, 1.0 / (IntAsDouble(order) - 1.0)));

That said, proposed change is indeed consistent with the paper.

[cgranade]

@guanghaolow: Thanks for the clarification; I think it would make sense given that to either add a remark to the API documentation comments for DecomposeIntoTimeStepsCA explaining the notation difference, or to make the change you suggested. In either case, I'll make an issue to track that improvement so that we can get this fix in in the meantime.