Я нашел решение, используя разложение ldl
L,d,_ = scipy.linalg.ldl(C)
d = np.diag(d).copy()
inds = d >= d.max()*1e-8
d = d[inds]
d = np.sqrt(d)
d.shape = (-1,1)
Q = d * L.T[inds]
loss = cp.norm(cp.matmul(Q, x))
Для разложения ldl требуется scipy> = 1.1.