Fixes #607

Introduction

Consider the singular value decomposition (SVD) of a complex-valued matrix $\mathbf{X} \in \mathbb{C}^{M \times N}$ truncated to rank $R$:
$$\mathbf{X} = \mathbf{U}\mathbf{S}\mathbf{V}^H, \quad \mathbf{U}\in\mathbb{C}^{M×R}, \mathbf{S}\in \mathbb{R}^{R×R}, \mathbf{V}^H \in \mathbb{C}^{R×N}$$
It is known that SVD is invariant under multiplication of each pair of singular vectors by a complex exponential $C_r = e^{j\phi(r)}$. In other words, for the $r$-th singular component we can replace: $\mathbf{u}_r \leftarrow C_r \cdot \mathbf{u}_r$ and $\mathbf{v}_r^H \leftarrow C_r^* \cdot \mathbf{v}_r^H$ without affecting the reconstruction accuracy.

Problem description

According to the documentation of the svd_flip function, its purpose is to eliminate the arbitrary phase of singular vectors by aligning them with the element of largest magnitude. In other words, the correction factor for each $r$ is chosen as:
$$C_r = \mathrm{sign}(z)^*,\quad z = \mathbf{u}_r(\mathrm{argmax}(|\mathbf{u}_r|))$$
Where $\mathrm{sign}(z) = z/|z|$.
But in current version the conjugate operation is missed.

Why it worked for numpy<2.0

Before Numpy 2.0, the definition of the complex sign function was different: $sign(z) = z/ \sqrt{zz}$ equvalent to (sign(x.real) + 0j if x.real != 0 else sign(x.imag) + 0j) which always evaluates to either $1$ or $-1$. Since both the left and right singular vectors were multiplied by the same $\pm 1$ factor, the overall decomposition remained consistent, and the missing conjugation had no visible effect.