According to Wikipedia:
A common convention is to list the singular values in descending order. In this case, the diagonal matrix $\Sigma$ is uniquely determined by $M$ (though the matrices $U$ and $V$ are not).
My question is, are $U$ and $V$ uniquely determined up to some equivalence relation (and what equivalence relation)?
Let $A = U_1 \Sigma V_1^* = U_2 \Sigma V_2^*$. Let us assume that $\Sigma$ has distinct diagonal elements and that $A$ is tall. Then
$$A^* A = V_1 \Sigma^* \Sigma V_1^* = V_2 \Sigma^* \Sigma V_2^*.$$
From this, we get
$$\Sigma^* \Sigma V_1^* V_2 = V_1^* V_2 \Sigma^* \Sigma.$$
Notice that $\Sigma^* \Sigma$ is diagonal with all different diagonal elements (that's why we needed $A$ to be tall) and $V_1^* V_2$ is unitary. Defining $V := V_1^* V_2$ and $D := \Sigma^* \Sigma$, we have
$$D V = V D.$$
Now, since $V$ and $D$ commute, they have the same eigenvectors. But, $D$ is a diagonal matrix with distinct diagonal elements (i.e., distinct eigenvalues), so it's eigenvectors are the elements of the canon basis. That means that $V$ is diagonal too, which means that
$$V = \operatorname{diag}(e^{{\rm i}\varphi_1}, e^{{\rm i}\varphi_2}, \dots, e^{{\rm i}\varphi_n}),$$
for some $\varphi_i$, $i=1,\dots,n$.
In other words, $V_2 = V_1 V$. Plug that back in the formula for $A$ and you get
$$A = U_1 \Sigma V_1^* = U_2 \Sigma V_2^* = U_2 \Sigma V^* V_1^* = U_2 V^* \Sigma V_1^*.$$
So, $U_2 = U_1 V$ if $\Sigma$ (and, in extension, $A$) is square nonsingular. Other options, somewhat similar to this, are possible if $\Sigma$ has zeroes on the diagonal and/or is rectangular.
If $\Sigma$ has repeating diagonal elements, much more can be done to change $U$ and $V$ (for example, one or both can permute corresponding columns).
If $A$ is not thin, but wide, you can do the same thing by starting with $AA^*$.
So, to answer your question: for a square, nonsingular $A$, there is a nice relation between different pairs of $U$ and $V$ (multiplication by a unitary diagonal matrix, applied in the same way to the both of them). Otherwise, you get quite a bit more freedom, which I believe is hard to formalize.
I'm going to provide a full characterisation of the set of SVDs for a given matrix $A$, using two different (but of course ultimately equivalent) kinds of formalisms. First a standard matrix formalism, and then using dyadic notation.
The TL;DR is that if $$A=UDV^\dagger=\tilde U \tilde D\tilde V^\dagger$$ for $U,\tilde U,V,\tilde V$ isometries and $D,\tilde D>0$ square strictly positive diagonal matrices, then we can safely assume $D=\tilde D$ by trivially rearranging the basis for the underlying space, and $\tilde V=V W$ and $\tilde U=U W$, with $W$ a block-diagonal unitary matrix such that $WD W^\dagger = D$. In particular, this means that $W$ can only mix subspaces that correspond to equal singular values. This characterises all and only the possible SVDs. So given any SVD, the freedom in the choice of other SVDs corresponds to the freedom in choosing these unitaries $W$ (how much freedom is that, depends in turn on the degeneracy of the singular values of $A$).
Unitaries vs isometries in the SVD — Note that I'm using here a convention where $D$ is square, and $V$ are isometries. This might appear slightly different than the other common method of taking $D$ generally rectangular and $U$ unitary, but it's essentially equivalent. If you take $D$ rectangular, and allow zeros on its diagonal, the corresponding columns of the unitary $U$ are irrelevant, and therefore we might just as well remove them, making $U$ an isometry and $D$ square with strictly positive diagonal. For a toy example of what I mean, consider the following SVD, written in the standard notation with unitaries and rectangular $D$: $$\begin{pmatrix}1 & 0 \\ 0 & 1 \\ 1 & 1\end{pmatrix} = \underbrace{\frac1{\sqrt6}\begin{pmatrix}1 & -\sqrt3 & -\sqrt2 \\ 1 & \sqrt3 & -\sqrt2 \\ 2 & 0 & \sqrt2\end{pmatrix}}_{U} \underbrace{\begin{pmatrix}\sqrt3 & 0\\0 & 1 \\ 0&0\end{pmatrix}}_{D} \underbrace{\frac{1}{\sqrt2}\begin{pmatrix}1&1\\-1&1\end{pmatrix}}_{V^\dagger}. $$ Clearly, the third column of $U$, which corresponds to the zero third row of $D$, is inconsequential. We can therefore more expressively write the SVD as $$\begin{pmatrix}1 & 0 \\ 0 & 1 \\ 1 & 1\end{pmatrix} = \underbrace{\frac1{\sqrt6}\begin{pmatrix}1 & -\sqrt3 \\ 1 & \sqrt3 \\ 2 & 0\end{pmatrix}}_{U} \underbrace{\begin{pmatrix}\sqrt3 & 0\\0 & 1\end{pmatrix}}_{D} \underbrace{\frac{1}{\sqrt2}\begin{pmatrix}1&1\\-1&1\end{pmatrix}}_{V^\dagger}, $$ by simply allowing $V$ to be an isometry rather than a unitary, and enforcing a strictly positive square diagonal $D$.
Consider the SVD of a given matrix $A$, in the form $A=UDV^\dagger$ with $D>0$ a square diagonal matrix with strictly positive entries, and $U,V$ isometries.
The question is, if you have $$A = UDV^\dagger = \tilde U \tilde D \tilde V^\dagger,$$ with $U,\tilde U,V,\tilde V$ isometries, and $D,\tilde D>0$ diagonal squared matrices, what does this imply for $\tilde U,\tilde D,\tilde V$? More specifically, can we somehow find an explicit relationship between them?
The first easy observation is that you must have $D=\tilde D$, modulo permutations of equal singular values (i.e. trivial rearrangements of the basis elements). This follows from $AA^\dagger=UD^2 U^\dagger = \tilde U \tilde D^2 \tilde U^\dagger$, which means that $D^2$ and $\tilde D^2$ both contain the eigenvalues of $AA^\dagger$. Because the spectrum is determined algebraically from $AA^\dagger$, the set of eigenvalues must be identical, and the singular values are by definition positive reals, we must have $D=\tilde D$.
The above reduces the question to: if $UDV^\dagger=\tilde U D\tilde V^\dagger$, with $U,\tilde U,V,\tilde V$ isometries, what can we say about $\tilde U,\tilde V$? To this end, we observe that the freedom in the choice of $U,V$ amounts to the different possible ways to decompose the subspaces associated to each distinct singular value.
More precisely, consider the subspace $V^{(d)}\equiv \{v: \, \|Av\|=d\}$ corresponding to a singular value $d$. We can then uniquely decompose the matrix as $A=\sum_d A_d$ where $A_d\equiv A \Pi_d$ and $\Pi_d$ is the projection onto $V^{(d)}$, and the sum is over all nonzero singular values $d$ of $A$. We can now observe that any and all SVDs of $A$ correspond to a choice of orthonormal basis for each $V^{(d)}$. Namely, for any such basis $\{\mathbf v_k\}$ we associate the partial isometry $V_d\equiv \sum_k \mathbf v_k \mathbf e_k^\dagger$. The corresponding basis for the image of $A_d$ is then determined as $\mathbf u_k= A \mathbf v_k$, and we then define the partial isometry $U_d\equiv \sum_k \mathbf u_k \mathbf e_k^\dagger$. Here, $\mathbf e_k$ denotes an orthonormal basis spanning the elements of $D$ corresponding to the singular value $d$. This procedure provides a decomposition $A_d= U_d D V_d^\dagger$, and therefore an SVD for $A$ itself by summing these. Any SVD can be constructed this way.
In conclusion, the freedom in choosing an SVD is entirely in the choice of bases $\{\mathbf v_k\}$ above. In matrix language, this means that given any SVD $A=UDV^\dagger$, all other SVDs can be written as $A=UW D (VW)^\dagger$ for some unitary $W$ that commutes with $D$. This commutation property is a concise way to state that $W$ is only allowed to mix vectors corresponding to the same singular value, that is, to the same eigenspace of $D$.
Let $$H \equiv \begin{pmatrix}1&1\\1&-1\end{pmatrix}.$$ This is a somewhat trivial example because $H$ is Hermitian and unitary. A standard SVD reads $$ H = \underbrace{\frac{1}{\sqrt2}\begin{pmatrix}1&1\\-1&1\end{pmatrix} }_{\equiv U} \underbrace{\begin{pmatrix}\sqrt2 & 0\\0&\sqrt2\end{pmatrix}}_{\equiv D} \underbrace{\begin{pmatrix}0&1\\1&0\end{pmatrix}}_{\equiv V^\dagger}.$$ In this case, we have two identical singular values. According to our discussion above, this means that we can apply a(ny) unitary transformation to the columns of $V$ and still obtain another SVD. That is, given any unitary $W$, $\tilde V\equiv V W$ gives another SVD for $H$. In this simple case, you can also observe this directly, as $D=\sqrt2 I$, and therefore $$H = UDV^\dagger= UD W \tilde V^\dagger = (UW) D \tilde V^\dagger,$$ hence $\tilde U\equiv UW$, $\tilde V\equiv VW$ give the alternative SVD $H=\tilde U D\tilde V^\dagger$, and all SVDs have this form.
For example taking $W=V$, we get $\tilde V=I$, and thus the decomposition $$H = \underbrace{\frac{1}{\sqrt2}\begin{pmatrix}1&1\\1&-1\end{pmatrix} }_{\equiv \tilde U} \underbrace{\begin{pmatrix}\sqrt2 & 0\\0&\sqrt2\end{pmatrix}}_{\equiv D} \underbrace{\begin{pmatrix}1&0\\0&1\end{pmatrix}}_{\equiv \tilde V^\dagger}.$$
Consider a simple non-squared case. Let $$A \equiv \begin{pmatrix}1&1\\1 & \omega \\ 1 & \omega^2\end{pmatrix}, \qquad \omega\equiv e^{2\pi i/3}.$$ This is again almost trivial because $A$ is an isometry, up to a constant. Still, we can write its SVD as $$A = \underbrace{\frac{1}{\sqrt3}\begin{pmatrix}1&1\\1&\omega\\1&\omega^2\end{pmatrix}}_{\equiv U} \underbrace{\begin{pmatrix}\sqrt3 & 0 \\0&\sqrt3\end{pmatrix}}_{\equiv D} \underbrace{\begin{pmatrix}0&1\\1&0\end{pmatrix}}_{\equiv V^\dagger}. $$ Notice that now $U,V$ are isometries but $U$ is not unitary, and that $D>0$ is squared. Per our results above, any SVD will have the form $A=\tilde U D \tilde V^\dagger$ with $\tilde V=V W$ for some unitary $W$.
for example, taking $W=V$ (we can do this, because here $V$ is also unitary), we get the alternative SVD $A=\tilde U D \tilde V^\dagger$ with $\tilde V=VW=VV=I$ and $\tilde U= U W^\dagger=UV=\frac1{\sqrt3}\begin{pmatrix}1&\omega&\omega^2\\1&1&1\end{pmatrix}^T$, that is, $$A = \underbrace{\frac{1}{\sqrt3}\begin{pmatrix}1&1\\\omega&1\\\omega^2&1\end{pmatrix}}_{\equiv \tilde U} \underbrace{\begin{pmatrix}\sqrt3 & 0 \\0&\sqrt3\end{pmatrix}}_{\equiv D} \underbrace{\begin{pmatrix}1&0\\0&1\end{pmatrix}}_{\equiv \tilde V^\dagger}.$$
Let's do an example with non-degenerate singular values. Let $$A = \begin{pmatrix}1& 2 \\ 1 & 2\omega \\ 1 & 2\omega^2\end{pmatrix}, \qquad \omega\equiv e^{2\pi i/3}.$$ This time the singular values are $\sqrt3$ and $2\sqrt3$. One SVD is easily derived as $$A = \underbrace{\frac{1}{\sqrt3}\begin{pmatrix}1&1\\1&\omega\\1&\omega^2\end{pmatrix}}_{\equiv U} \underbrace{\begin{pmatrix}\sqrt3 & 0 \\0&2\sqrt3\end{pmatrix}}_{\equiv D} \underbrace{\begin{pmatrix}1&0\\0&1\end{pmatrix}}_{\equiv V^\dagger}.$$ However, in this case there is much less freedom in choosing other SVDs, because these must correspond to $\tilde V=VW$ where $W$ only mixes columns of $V$ corresponding to the same values of $D$. In this case $D$ is non-degenerate, thus $W$ must be diagonal with unimodular entries, and therefore the full set of SVDs must correspond to $W=\begin{pmatrix}e^{i\alpha}&0\\0&e^{i\beta}\end{pmatrix}$, that is, $$A = \underbrace{\frac{1}{\sqrt3}\begin{pmatrix}e^{-i\alpha}&e^{-i\beta}\\e^{-i\alpha}&\omega e^{-i\beta}\\e^{-i\alpha}&\omega^2 e^{-i\beta}\end{pmatrix}}_{\equiv \tilde U} \underbrace{\begin{pmatrix}\sqrt3 & 0 \\0&2\sqrt3\end{pmatrix}}_{\equiv D} \underbrace{\begin{pmatrix}e^{i\alpha}&0\\0&e^{i\beta}\end{pmatrix}}_{\equiv \tilde V^\dagger}.$$ All SVDs will look like this, for some $\alpha,\beta\in\mathbb{R}$. Permuting the elements of $D$ we can obtain SVDs which look superficially different but are ultimately equivalent to the above.
The SVD of an arbitrary matrix $A$ can be written in dyadic notation as $$A=\sum_k s_k u_k v_k^*,\tag A$$ where $s_k\ge0$ are the singular values, and $\{u_k\}_k$ and $\{v_k\}_k$ are orthonormal sets of vectors spanning $\mathrm{im}(A)$ and $\ker(A)^\perp$, respectively. The connection between this and the more standard way of writing the SVD of $A$ as $A=UDV^\dagger$ is that $u_k$ is the $k$-th column of $U$, and $v_k$ is the $k$-th column of $V$.
If $A$ is nondegenerate, the only freedom in the choice of vectors $u_k,v_k$ is their global phase: replacing $u_k\mapsto e^{i\phi}u_k$ and $v_k\mapsto e^{i\phi}v_k$ does not affect $A$.
On the other hand, when there are repeated singular values, there is additional freedom in the choice of $u_k,v_k$, similarly to how there is more freedom in the choice of eigenvectors corresponding to degenerate eigenvalues. More precisely, note that (A) implies $$AA^\dagger=\sum_k s_k^2 \underbrace{u_k u_k^*}_{\equiv\mathbb P_{u_k}}, \qquad A^\dagger A=\sum_k s_k^2 \mathbb P_{v_k}.$$ This tells us that, whenever there are degenerate singular values, the corresponding set of principal components is defined up to a unitary rotation in the corresponding degenerate eigenspace. In other words, the set of vectors $\{u_k\}$ in (A) can be chosen as any orthonormal basis of the eigenspace $\ker(AA^\dagger-s_k^2)$, and similarly $\{v_k\}_k$ can be any basis of $\ker(A^\dagger A-s_k^2)$.
However, note that a choice of $\{v_k\}_k$ determines $\{u_k\}$, and vice-versa (otherwise $A$ wouldn't be well-defined, or injective outside its kernel).
A choice of $U$ uniquely determines $V$, so we can restrict ourselves to reason about the freedom in the choice of $U$. There are twe main sources of redundancy:
Finally, we should note that the former point is included in the latter, which therefore encodes all of the freedom allowed in choosing the vectors $\{v_k\}$. This is because multiplying the elements of an orthonormal basis by phases does not affect its being an orthonormal basis.
I will complete the proof of @Vedran for the case when there exist repeating eigenvalues, which would justify what @glS have said. Let $A = U \Sigma V^T = U^{'} \Sigma {V^{'}}^T$ be a matrix with real entries - the case with complex entries is similar. Then, $$A^T A = V \Sigma^T \Sigma V^{T} = V^{'} \Sigma^T \Sigma {V^{'}}^T.$$
From this, we get $$\Sigma^T \Sigma V^T V^{'} = V^T V^{'} \Sigma^T \Sigma.$$ Defining the square matrix $Q$ as $Q = V^T V^{'}$, we have $Q^T Q = (V^T V^{'})^T V^T V^{'} = I$, and similarly, $Q Q^T = I.$ Hence, $Q$ is an orthogonal matrix that satisfies the Sylvester equation $$Q\Sigma^T\Sigma - \Sigma^T \Sigma Q = 0.\tag{1}$$
Aiming to simplify the Sylvester equation (1) a little, counting the multiplicities, we can write $\Sigma^T \Sigma $ $= \sigma_1^2 I_{n_1} \oplus \sigma_2^2 I_{n_2} \oplus \cdots $ $ \oplus \sigma_{k}^2 I_{n_k} $ $= \text{diag}(\sigma_1^2 I_{n_1}, \sigma_2^2 I_{n_2},$ $ \cdots, \sigma_{k}^2 I_{n_k}),$ with $\sigma_i=\sigma_j$ iff $i=j$ and $n_i$ the multiplicity of $\sigma_{i}$.
Now, writing $Q$ in blocks conformally to $\Sigma^T \Sigma$, i.e., with $Q \Sigma^T \Sigma$ and $\Sigma^T \Sigma Q$ making sense, we have a new system of Sylvester equations $$\sigma_i^2 Q_{ij} I_{n_i} - \sigma_j^2 I_{n_j} Q_{ij}=0,$$ for each $1\leq i,j\leq k$. This means that, since both matrices $\sigma_i^2 I_{n_i}$ and $\sigma_j^{2}I_{n_j}$ do not share any of its eigenvalues for $i\not=j,$ by the result discussed here, if $i\not= j$, we have necessarily that $Q_{ij} = 0$. This means that $V^T V^{'} = \text{diag} (Q_{11}, Q_{22},\cdots, Q_{kk})$ for orthogonal matrices $Q_{ii}$, meaning that $$V' = V \text{diag}(Q_{11}, Q_{22},\cdots, Q_{kk})=V Q.\tag{2}$$
We can repeat the argument above for the matrix $A A^T$ and conclude that there exist an orthogonal matrix $\bar{Q}$ such that $$\bar{Q} = \text{diag}(\bar{Q}_{11}, \bar{Q}_{22},\cdots, \bar{Q}_{kk})\tag{3}$$ and $U' = U \bar{Q},$ with the matrices $\bar{Q}_{ii}$ of the same size of $Q_{ii}$ whenever $\sigma_i\not=0$ - this is possible because the matrices $A^T A$ and $A A^T$ have the same eigenvalues with the same multiplicity, except possibly the null eigenvalue. Taking into account that $A = U \Sigma V^T = U^{'} \Sigma {V^{'}}^T$, we would be able to conclude $\Sigma = U^T U^{'} \Sigma {V^{'}}^T V = U^T U \bar{Q} \Sigma Q^T V^T V,$ which simplifies to $$\Sigma = \bar{Q} \Sigma Q^T. $$ Lastly, considering the nonzero blocks of $\Sigma$ associated with the matrices $Q_{ii}$ and $\bar{Q}_{ii}$, we are able to conclude that $\bar{Q}_{ii} = Q_{ii}$ in the expression (3), whenever $\sigma_i\not=0$.
$\newcommand\R{\mathbb{R}}$I see that this question has been asked a few times. I'd like to provide a simple answer to this.
The singular values $s_1, \dots, s_k$ of an $n$-by-$m$ matrix $M$ are the square roots of the positive eigenvalues of $M^*M$. Let $d_1, \dots, d_k$ be the respective dimensions of the eigenspaces, and $d_{k+1}$ be the dimension of the kernel. Observe that $r = d_1+\cdots + d_k$ is the rank of $M$.
If $$ M = U\Sigma V^* $$ is a singular decomposition, then $M^*M = V\Sigma^2V^*$, and therefore the columns of $V$ form a unitary basis of eigenvectors $(v_1, \dots, v_m)$. $V$ is clearly unique up to unitary transformations of each eigenspace.
Now, for each $1 \le j \le k$, let $w_j = Mv_j$ and observe that for any $1 \le i,j \le r$, $$ \langle w_i,w_j\rangle = \langle Mv_i,Mv_j\rangle = \langle v_i,M^*Mv_j\rangle = \lambda_j\langle v_i,v_j\rangle = \lambda_j\delta_{ij}, $$ where the $\lambda_j$ are the singular values listed with the appropriate multiplicities. Therefore, $$ u_j = \frac{w_j}{\lambda_j},\ 1 \le j \le r $$ are a unitary basis of the image of $M$. They can be extended to a unitary basis $(u_1, \dots, u_n)$ of $\R^n$. It is now easy to show that the matrix $U$ whose columns are $u_1, \dots, u_n$ satisfies the equation $$ M = U\Sigma V^*.$$ It is also easy to show that conversely, if $$ M = U\Sigma V^*, $$ then, for each $1 \le j \le r$, then $$ \lambda_j u_j = Mv_j,$$ and $u_{r+1}, \dots, u_n$ are a unitary basis of the orthogonal complement of the image of $M$. It follows that given any $V$ satisfying the conditions above, the first $r$ columns of $U$ are uniquely determined and the last $n-r$ columns can be any unitary basis of $(\operatorname{image} M)^\perp$.
From all this, we see that $U$ and $V$ are unique up to unitary transformations of the eigenspaces of $M^*M$ and unitary transformations of $(\operatorname{image} M)^\perp$.
