Methods for estimating a 2D transformation without knowing the point matchings $\mathbf{p}=[x_P,y_P]^T$ and $\mathbf{q}$ are 2D points.
The transformation $f:\mathbb{R}^2\rightarrow\mathbb{R}^2$ maps $\mathbf{p}$ to $\mathbf{q}$ and has a vector of parameters $\boldsymbol{\theta}$:
$\mathbf{q}=f(\mathbf{p},\boldsymbol{\theta})$
For example, in a translation, $\boldsymbol{\theta}=[x_T,y_T]$ and $f$ is:
$x_Q=x_P+x_T$
$y_Q=y_P+y_T$
Now I have a list $S$ of 2D points $S=(\mathbf{p}_1,\mathbf{p}_2,\cdots,\mathbf{p}_N)$ and a list $D=(\mathbf{q}_1,\mathbf{q}_2,\cdots,\mathbf{q}_N)$,
$\mathbf{q}_1=f(\mathbf{p}_1,\boldsymbol{\theta})$
$\mathbf{q}_2=f(\mathbf{p}_2,\boldsymbol{\theta})$
$\vdots$
$\mathbf{q}_N=f(\mathbf{p}_N,\boldsymbol{\theta})$
If I know that the first point in $S$ is mapped to the first point in $D$, the second is mapped to the second, and so on, and assuming to know that $f$ is a translation/rotation/mirroring/etc., I will for example try to estimate $\boldsymbol{\theta}$ in this way:
$\underset{\boldsymbol\theta}{\operatorname{argmin}} \displaystyle\sum\limits_{i=1}^N (\mathbf{q}_i-f(\mathbf{p}_i,\boldsymbol{\theta}))^2$
Then someone gives me a random permutation $S'$ and $D'$ of $S$ and $D$ and does not tell me the permutation: now I can no longer apply the above minimization because I do not know how to match the $\mathbf{p}$s in $S'$ with the $\mathbf{q}$s in $D'$.
Is there any method/algorithm which simultaneously gives me the matchings and the $\boldsymbol{\theta}$?
I am interested in the general case even if my case study is the Reflection
 where $\mathbf{\theta}$ is given by $\mathbf{x}_0$ and $\mathbf{n}$ and $f$ is
$\mathbf{q}=-\mathbf{p}+2\:\mathbf{x}_0+2\:\hat{\mathbf{n}}[(\mathbf{x}_1-\mathbf{x}_0)\cdot\hat{\mathbf{n}}]$
 A: For the reflection case.
Given $P=\{p_k\},\ k=1,\cdots ,n$ a set of points containing their mirror image, we have that $P$ has at least one symmetry axis which characterizes the mirroring line. The procedure to determine the mirroring line is linked to the eigenvectors determination for the matrix
$$
I = \left(\begin{array}{cc}
\sum_{k=1}^n (x_k-x_g)^2 & -2\sum_{k=1}^n (x_k-x_g)(y_k-y_g)\\
-2\sum_{k=1}^n (x_k-x_g)(y_k-y_g) & \sum_{k=1}^n (y_k-y_g)^2
\end{array}\right)
$$
where $p_k = (x_k, y_k)$ and $g = (x_g,y_g) = \frac 1n\sum_{k=1}^n p_k$.
NOTE
One of the eigenvectors can be used as direction for the mirroring line and $g$ is contained into this line, thus characterizing the reflection. $P$ can have many symmetry axis.
A: The problem you are describing is called point cloud registration. The objective of point cloud registration is to find the transformation between point clouds.
Using your notation, we have a set of point $P = \{ \mathbf{p}_0, \mathbf{p}_1, \cdots, \mathbf{p}_N \}$ and a set of points $Q = \{ \mathbf{q}_0, \mathbf{q}_1, \cdots, \mathbf{q}_N \}$. This problem is typically written as follows:
$$
\begin{align}
R, t = \underset{R, t}{\operatorname{argmin}} \sum_{i=0}^{N-1}|| \mathbf{p}_i - R \mathbf{q}_i - t ||^{2}
\end{align}
$$
where $R$ is a rotation matrix, $t$ is a translation vector, and $\mathbf{p}_i$ and $\mathbf{q}_i$ are point correspondences with $N$ points. This can also be written using homogeneous transformations, but the problems are equivalent. If the point correspondences are known, then the optimal transformation can be obtained using singular value decomposition, but typically, these are unknown (as you mentioned).
If unknown, the problem can be solved using the Iterative Closest Point (ICP) algorithm, which is probably the most commonly used algorithm for this task. However, several algorithm exist that usually perform better than ICP such as TEASER, which is a certifiable algorithm for point cloud registration.
Note, I believe the mentioned algorithms do not consider reflections. The problem as I described solves for $R$ and $t$ where $R$ is a member of the special orthogonal group, which is the group of all orthogonal matrices of determinant $+1$ (i.e., group of all rotation matrices). Some sorts of transformation have determinant $-1$, which represent reflections, but the problem described above does nt consider determinant $-1$.
