Hadamard transform efficient

[AGaliciaMartinez]

Description
Improved efficiency of the hadamard gate creation time as it was using python for loops which are terribly slow. For N=10 it goes down from 522ms to roughly 7ms.

Changelog
Improved hadamard_transform efficiency (70x faster for N=10).

[quantshah]

@BoxiLi How should we handle such changes to qutip/qip in v4.6 since for qutip 5 this will be deprecated. Should these improvements be made directly to qutip-qip? We do not want to end up with three different versions (qutip 4.x, qutip 5 and qutip-qip).

[hodgestar]

Left one tiny suggestion, but this looks great -- easier to read and 70x faster. Thank you for adding the test.

What made you trip over this?

[hodgestar]

@quantshah @BoxiLi I assume the plan is to backport this to qutip-qip (as with the other qutip.qip changes in 4.X)?

@AGaliciaMartinez You're welcome to open a copy of this change in the qutip-qip repo once it's merged if you like? That would save @BoxiLi doing it (and it's good to have more eyes on qutip-qip).

[AGaliciaMartinez]

What made you trip over this?

I was using this function with N=13 for some testing code I have and realised it took too long to create the matrix.

[coveralls]

Coverage decreased (-0.01%) to 66.807% when pulling 311ba24 on AGaliciaMartinez:hadamard_transform_efficient into 5060a7d on qutip:master.

[BoxiLi]

As it is quite a short (but efficient!) change, we can keep it both in qutip and in qutip-qip.

Also, I remember this function is also used in the control module (also for testing I think). So we need it in qutip anyway.

Eventually, I think it is best to keep these commonly used function in qutip.operator and import it to qutip-qip

[jakelishman]

There's another interesting technique that's sometimes applicable in these sort of repeating cases. Rather than passing a list of the same elements to tensor, you can also consider building up the tensor product manually, going up in powers of two. So you go along with something like:

out, tmp = H, H
for i in n.bit_length():
    tmp = tensor(tmp, tmp)
    if n & (1 << i):
        out = tensor(out, tmp)

(very very approximately - I've almost certainly got the indexing wrong). The idea is that you only perform lg(n) tensor products.

Now, this might not actually give you much of a speed up in this case, because the calculation is still going to be dominated by the final tensor product, since the output matrix keeps getting larger. But in things like integer matrix powers, bigint calculations, or binomial expansions, this type of thing can be super useful! QuTiP 5 uses an algorithm like this for matrix powers. (I imagine SciPy does too, but for all the normal reasons, we don't use their implementation.)

edit: here it is:

qutip/qutip/core/data/pow.pyx

Lines 17 to 39 in 559387e

	cpdef CSR pow_csr(CSR matrix, unsigned long long n):
	if matrix.shape[0] != matrix.shape[1]:
	raise ValueError("matrix power only works with square matrices")
	if n == 0:
	return csr.identity(matrix.shape[0])
	if n == 1:
	return matrix.copy()
	# We do the matrix power in terms of powers of two, so we can do it
	# ceil(lg(n)) + bits(n) - 1 matrix mulitplications, where `bits` is the
	# number of set bits in the input.
	#
	# We don't have to do matrix.copy() or pow.copy() here, because we've
	# guaranteed that we won't be returning without at least one matrix
	# multiplcation, which will allocate a new matrix.
	cdef CSR pow = matrix
	cdef CSR out = pow if n & 1 else None
	n >>= 1
	while n:
	pow = matmul_csr(pow, pow)
	if n & 1:
	out = pow if out is None else matmul_csr(out, pow)
	n >>= 1
	return out