ParamResolver is slow

[peterse]

Reconstructing a circuit op-wise and replacing all parametrized ops with new ops containing parameter updates is 10-50x faster than using cirq.ParamResolver with sympy symbols...

Example code to reproduce:

from timeit import default_timer as timer
import numpy as np
import sympy
import cirq


def update_params(circuit, params):
    """Competitor method to cirq.ParamResolver."""
    new_op_tree = []
    for op, param in zip(circuit.all_operations(), params):
        new_op_tree.append(op.gate._with_exponent(param/np.pi)(*op.qubits))
    return cirq.Circuit.from_ops(new_op_tree)


trials = 100
for depth in [10, 15, 20]:
    sympy_circuit = cirq.Circuit.from_ops([cirq.Rx(sympy.Symbol(str(k)))(cirq.LineQubit(0)) for k in range(depth)])
    random_params = np.random.randn(trials, depth)
    # time twenty runs
    start = timer()
    for j in range(trials):
        resolver = dict(zip([str(k) for k in range(depth)], random_params[j]))
        wf1 = cirq.Simulator().simulate(sympy_circuit, param_resolver=resolver).final_state
    end = timer() - start
    print(f"{depth} parameters, {trials} trials using Sympy+ParamResolver: {end} seconds")

    start = timer()
    for j in range(trials):
        float_circuit = update_params(sympy_circuit, random_params[j])
        wf2 = cirq.Simulator().simulate(float_circuit).final_state
    end = timer() - start
    print(f"{depth} parameters, {trials} trials using reconstructed circuit: {end} seconds")

produced (cirq v0.5.0, Windows 10 + core i7 gen 7 processor)

>>> 10 parameters, 100 trials using Sympy+ParamResolver: 2.408036000095308 seconds
>>> 10 parameters, 100 trials using reconstructed circuit: 0.1671589999459684 seconds
>>> 15 parameters, 100 trials using Sympy+ParamResolver: 4.347879000008106 seconds
>>> 15 parameters, 100 trials using reconstructed circuit: 0.25207799999043345 seconds
>>> 20 parameters, 100 trials using Sympy+ParamResolver: 7.1194350000005215 seconds
>>> 20 parameters, 100 trials using reconstructed circuit: 0.31734399998094887 seconds

Some plots of how this scales with large numbers of parameters (generated from a different script using the same update_params method):

[dabacon]

Interesting. We should narrow this down to whether it is something in the Simulator or the ParamResolver or the protocol. A first glance at code suggest that hashing in ParamResolver, that doesn't seem like a great idea.

[Ashalynd]

Some thoughts:

1. ParamResolver initialization might be potentially heavy, but afaics it is not supposed to be re-initialized in the recursive calls.
2. Recursion itself might potentially be a problem, which is a known issue in Python. The winning solution doesn't do recursion for parameter resolving.

[peterse]

@Ashalynd could you say a little more about how recursion is used in the resolver?

[alexandrupaler]

The following line where sympy.subs

Cirq/cirq/study/resolver.py

Line 74 in 57b0beb

v = value.subs(self.param_dict)

consumes most of the time

sympy.subs is known to be very slow, and is called very often (and includes, as @Ashalynd said, a lot of recursions).

After having a look in subs I came up wit the following code to replace the line mentioned above

v = value
for pr in value.args:
    str_pr = str(pr)
    if str_pr in self.param_dict:
        v = v.subs(str_pr, self.param_dict[str_pr])

On my computer it is now 3 times faster, but still approx. 30 times slower than the method without sympy.

Would a solution be to use lambdification?

[dabacon]

Here is the page on sub being slow and lambification

[Strilanc]

We could do this:

def Rz(rads: value.TParamVal) -> ZPowGate:
    """Returns a gate with the matrix e^{-i Z rads / 2}."""
    if isinstance(rads, sympy.Symbol):
        return ZPowGate(exponent=FastOverPi(rads), global_shift=-0.5)
    pi = sympy.pi if protocols.is_parameterized(rads) else np.pi
    return ZPowGate(exponent=rads / pi, global_shift=-0.5)

class FastOverPi(sympy.Mul):
    def __init__(self, v: sympy.Symbol):
        super().__init__(v, 1 / sympy.pi)
        self.v = v

    def subs(self, *args):
        if isinstance(args[0], dict):
            return args[0].get(str(self.v), self.v) / np.pi
        return super().subs(*args)

[peterse]

@karlunho

[alexandrupaler]

sympy is generally slow, and everything that touches sympy becomes slow.

Would it be possible to have simple parameterized gates (e.g. the Rz rotations uses sympy only to multiply a number by pi) without sympy? If the symbol in the constructor stands just for a number, sympy is overkill, if the symbol is an expression per se, maybe the sympy overhead is good price.

[dstrain115]

I took a stab at short-circuiting sympy for common cases for PR #2394 . It does not help the above code as written, but, if flattened or if using XPowGate (which doesn't add the "x/pi" formula), this speeds up the code significantly.

10 parameters, 100 trials using flattened circuit: 0.1600349259097129 seconds
10 parameters, 100 trials using xpow circuit: 0.17739228205755353 seconds
10 parameters, 100 trials using reconstructed circuit: 0.16374025191180408 seconds:
15 parameters, 100 trials using flattened circuit: 0.2257283239159733 seconds
15 parameters, 100 trials using xpow circuit: 0.2516175350174308 seconds
15 parameters, 100 trials using reconstructed circuit: 0.23388241999782622 seconds:
20 parameters, 100 trials using flattened circuit: 0.2913768459111452 seconds
20 parameters, 100 trials using xpow circuit: 0.32354101818054914 seconds
20 parameters, 100 trials using reconstructed circuit: 0.29962173686362803 seconds:

The changed code is below, for reference:

trials = 100
for depth in [10, 15, 20]:
    global sympy_circuit
    global resolver
    sympy_circuit1 = cirq.Circuit.from_ops([cirq.Rx(sympy.Symbol(str(k)))(cirq.LineQubit(0)) for k in range(depth)])
    sympy_circuit2 = cirq.Circuit.from_ops([cirq.XPowGate(exponent=sympy.Symbol(str(k)))(cirq.LineQubit(0)) for k in range(depth)])
    sympy_circuit, expr_map = cirq.flatten(sympy_circuit1)
    random_params = np.random.randn(trials, depth)
    resolver = expr_map.transform_params(dict(zip([str(k) for k in range(depth)], random_params[0])))
    # time twenty runs
    start = timer()
    for j in range(trials):
        wf = cirq.Simulator().simulate(sympy_circuit, param_resolver=resolver).final_state
    end = timer() - start
    print(f"{depth} parameters, {trials} trials using flattened circuit: {end} seconds")

    # time twenty runs
    start = timer()
    for j in range(trials):
        resolver = dict(zip([str(k) for k in range(depth)], random_params[j]))
        wf = cirq.Simulator().simulate(sympy_circuit2, param_resolver=resolver).final_state
    end = timer() - start
    print(f"{depth} parameters, {trials} trials using xpow circuit: {end} seconds")

    start = timer()
    for j in range(trials):
        float_circuit = update_params(sympy_circuit, random_params[j])
        wf2 = cirq.Simulator().simulate(float_circuit).final_state
    end = timer() - start
    print(f"{depth} parameters, {trials} trials using reconstructed circuit: {end} seconds:")

I'm still facing a problem. In my case, new theta values are generated in a different scope and delivered to the ParamSolver. The next value of thetas would be unknown to me. I will only know it after I ParamSolver and generate the unitary of the circuit for the current thetas. So for me, it has to be one theta set at a time. Any suggestions?

[dstrain115]

@aminjahanpour Can you provide a simplified example of the problem? Otherwise, I am not sure I can reproduce your issue.

One suggestion would be to reduce the use of formulas as much as possible, since sympy formula resolution is extremely slow. If you are using the Rz/Ry/Rx gates, you may want to consider using ZPowGate,YPowGate,XPowGate instead and divide by pi yourself rather than rely on sympy to resolve it.