Description
For my application, Qobj.ptrace() was a bottleneck, so I added a specialization for kets that is significantly faster.
The new implementation uses np.einsum to combine the computation of the outer product and the partial trace into one operation, instead of doing these two operations sequentially.

from qutip import rand_ket
from time import perf_counter

state = rand_ket([2, 3, 4, 5]*2, distribution="fill")
keep_dims = [1, 3, 7]

t0 = perf_counter()
reduced_dm2 = state.ptrace(keep_dims) # Uses new implementation
t1 = perf_counter()
print("New impementation:", t1 - t0) # ~ 0.01s on 8th-gen Intel quad-core


t0 = perf_counter()
reduced_dm1 = state.proj().ptrace(keep_dims) # Old implementation
t1 = perf_counter()
print("Old impementation:", t1 - t0) # ~ 6s on 8th-gen Intel quad-core

assert (reduced_dm1 - reduced_dm2).norm("tr") < 1e-10

In this example, the new implementation is around 600x faster on my machine.

The new implementation can also handle significantly larger states:
rand_ket([2, 3, 4, 5]*3, distribution="fill").ptrace(keep_dims) fails using the old implementation because the intermediate density matrix is too large, but the new implementation computes the result in less than 0.1s.

Note: I did not use qutip.einsum because this only seems to wrap the string-indexed version of np.einsum and not the array-indexed version. The latter is more useful is cases like the present one, since its easier do generate the einsum subscripts programatically.