Here is a more detailed description of what I wrote in the comments. Ideally such a description should come with a ton of graphics illustrating all of the discussed ideas. I did not have the time to produce said images, so my answer is lacking in that respect.
Summary For Implementation Purposes
Let the vertices of the first triangle be $v_0, v_1, v_2\in\mathbb{R}^3$ and the vertices of the second triangle be $w_0, w_1, w_2\in\mathbb{R}^3$. For a unique solution we require that the first triangle is non-degenerate: it must have non-zero area (we cannot have the 3 vertices lying on a single line or in a single point).
Compute $e_1 = v_1-v_0, \, e_2 = v_2-v_0$ and $f_1 = w_1-w_0, \, f_2 = w_2 - w_0$.
Option 1
Compute $e_3 = e_1 \times e_2$ and $f_3 = f_1 \times f_2$ where $\times$ is the cross product. Then form the $3\times 3$ matrices $E$ and $F$ using the $e_i$ and $f_i$ as column vectors:
$$
\begin{gather}
E = \begin{bmatrix} e_1 & e_2 & e_3 \end{bmatrix}, \quad F = \begin{bmatrix} f_1 & f_2 & f_3\end{bmatrix}.
\end{gather}
$$
Finally, compute $M = FE^{-1}$ where $E^{-1}$ is the inverse of $E$, and compute $m_3 = w_0 - Mv_0$. Set $M$ as the $3\times 3$ part of your $4\times 4$ matrix $M'$ and $m_3$ as its last column:
$$
\begin{equation}
M' = \begin{bmatrix} m_0 & m_1 & m_2 & m_3 \\ 0 & 0 & 0 & 1 \end{bmatrix}.
\end{equation}
$$
Option 2
Form the $3\times 2$ matrices $E$ and $F$:
$$
\begin{equation}
E = \begin{bmatrix} e_1 & e_2 \end{bmatrix}, \quad F = \begin{bmatrix} f_1 & f_2 \end{bmatrix}.
\end{equation}
$$
The form of the solution is similar to the one from Option 1:
$$
\begin{gather}
M = FE^{+}, \quad m_3 = w_0 - FE^{+}v_0,\quad E^{+} = (E^TE)^{-1}E^T \\
M' = \begin{bmatrix} m_0 & m_1 & m_2 & m_3 \\ 0 & 0 & 0 & 1 \end{bmatrix}.
\end{gather}
$$
Differences Between Option 1 and Option 2
The two approaches (Option 1 and Option 2) agree when applied to points $p$ that lie in the plane of the triangle $v_0, v_1, v_2$, however they differ if $p$ is not in this plane.
Option 1 rescales space perpendicular to the triangles $v_0, v_1, v_2$ and $w_0, w_1, w_2$ according to the ratio $r = \frac{area(w_0, w_1, w_2)}{area(v_0,v_1,v_2)}$. For example, if $r = 2$ (the second triangle has twice larger area than the first) and $p$ was $3$ units perpendicularly away from the first triangle's plane, then its transformed version would be $3\times r = 3\times 2 = 6$ units perpendicularly away from the second triangle.
If $p$ is not inside the plane in which $v_0, v_1, v_2$ lies, then Option 2 projects $p$ on this plane, and only then does it transform it. That is, you would always get a transformed point in the plane of the second triangle $w_0,w_1,w_2$ regardless of what point you feed in, since it would first get projected on the plane of the first triangle. For example if you had a point $p$ that is $3$ units perpendicularly away from the plane of the first triangle, then it will end up $0$ units perpendicularly away from the plane of the second triangle (i.e. it will be in it).
Introduction
The following sections contain a more detailed treatment of the problem along with guidelines on ideas that may be used to derive the final solution.
I'll start with a simplified case for clarity. It turns out that studying affine transformations between tetrahedra in 3-D is simpler than studying affine transformations between triangles in 3-D. Studying transformations between triangles is simplest in 2-D for instance, and for $n$-D it is simplest to study transformations between $n$-simplices (a $2$-simplex is a triangle, and a $3$-simplex is a tetrahedron). This is due to the fact that an $n$-dimensional vector space is uniquely represented through exactly $n$ basis vectors. Consequently, studying transformations between triangles (a 2-D object) in 3-D results in some difficulties compared to studying tetrahedra, since in 3-D the triangle case leaves out one of the basis vectors.
Tetrahedron Case: Formulation
Let $v_0,v_1,v_2, v_3 \in\mathbb{R}^3$ be the vertices of a non-degenerate tetrahedron $\mathcal{T}_{\mathcal{V}}$ in 3-D, and let $w_0, w_1, w_2, w_3\in\mathbb{R}^3$ be the vertices of another (not necessarily non-degenerate) tetrahedron $\mathcal{T}_{\mathcal{W}}$. We want an affine (linear + offset) function $g:\mathcal{T}_{\mathcal{V}}\to\mathcal{T}_{\mathcal{W}}$ that maps the first tetrahedron to the second one. As a special case the corresponding vertices must map to each other:
$$
\begin{gather}
g(v_0) = w_0, \quad
g(v_1) = w_1, \quad
g(v_2) = w_2, \quad
g(v_3) = w_3.
\end{gather}
$$
Since $g(p)$ is affine and acts on 3-D points, then it follows from the definition that it can be represented by a multiplication with a 3x3 matrix and a translation:
$$
\begin{equation}
g(p) = Mp + m_3, \quad M = \begin{bmatrix} m_0 & m_1 & m_2 \end{bmatrix}, \quad m_0, m_1, m_2, m_3 \in\mathbb{R}^3.
\end{equation}
$$
That is, we have a 3x3 matrix $M$ with 3-D column vectors $m_0, m_1, m_2$ and one translation/offset 3-D vector $m_3$. Rewriting the constraints using the above form of $g$ yields:
$$
\begin{align}
g(v_0) = w_0 &\iff m_0 v_{0,x} + m_1 v_{0,y} + m_2 v_{0,z} + m_3 = w_0 \\
g(v_1) = w_1 &\iff m_0 v_{1,x} + m_1 v_{1,y} + m_2 v_{1,z} + m_3 = w_1 \\
g(v_2) = w_2 &\iff m_0 v_{2,x} + m_1 v_{2,y} + m_2 v_{2,z} + m_3 = w_2 \\
g(v_3) = w_3 &\iff m_0 v_{3,x} + m_1 v_{3,y} + m_2 v_{3,z} + m_3 = w_3.
\end{align}
$$
Tetrahedron Case: Solution
The above is a system with $4\times 3=12$ unknowns (4 $m_i$ vectors each with 3 components) and $4\times 3=12$ equations. From linear algebra we know that as long as the equations are not linearly dependent this system has a unique solution. We can express $m_0$ from the first equation:
$$
\begin{equation}
m_3 = w_0 - (m_0 v_{0,x} + m_1 v_{0,y} + m_2 v_{0,z}).
\end{equation}
$$
Plugging this expression in the remaining equations reduces the system to 9 unknowns:
$$
\begin{align}
m_0 e_{1,x} + m_1 e_{1,y} + m_2 e_{1,z} &= f_1 \\
m_0 e_{2,x} + m_1 e_{2,y} + m_2 e_{2,z} &= f_2 \\
m_0 e_{3,x} + m_1 e_{3,y} + m_2 e_{3,z} &= f_3,
\end{align}
$$
where $e_1 = v_1-v_0, \, e_2 = v_2-v_0, \, e_3 = v_3-v_0$ and $f_1 = w_1-w_0, \, f_2 = w_2-w_0, \, f_3 = w_3-w_0$. This can be written compactly in matrix form:
$$
\begin{gather}
M E = F, \, E = \begin{bmatrix} e_1 & e_2 & e_3 \end{bmatrix}, \, F = \begin{bmatrix} f_1 & f_2 & f_3 \end{bmatrix}\\
\implies \\
M = E^{-1}F, \, m_3 = w_0 - Mv_0.
\end{gather}
$$
The matrix $E^{-1}$ exists provided that its columns are linearly independent. This is the case when the tetrahedron with vertices $v_0, v_1, v_2, v_3$ is non-degenerate: the points do not coincide, they do not lie on a single line, they do not lie in a single plane, i.e. it's a proper pyramid with non-zero volume (not a squished one).
We can form a $4\times 4$ matrix by using $M$ as the $3\times 3$ part, $m_3$ as the last column, and setting $\begin{bmatrix} 0 & 0 & 0 & 1 \end{bmatrix}$ in the last row. Then the transformation takes the form:
$$
\begin{equation}
f(p) = q \iff \begin{bmatrix} m_0 & m_1 & m_2 & m_3 \\ 0 & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} p \\ 1 \end{bmatrix} = \begin{bmatrix} q \\ 1 \end{bmatrix}
\end{equation}
$$
Note also that we could have taken any of the vertices as origin instead of $v_0, w_0$:
$$
\begin{gather}
m_3 = w_i - (m_0 v_{i,x} + m_1 v_{i,y} + m_2 v_{i,z}) \\
e_j = v_j - v_i, \, f_j = w_j - w_i, \, j\ne i,
\end{gather}
$$
so the choice of $v_0,w_0$ as origins in the two coordinate systems was arbitrary.
Triangle Case: Formulation
The setting of transforming between two triangles in 3-D is similar to the tetrahedra case, however $v_3$ and $w_3$ are not provided. We have:
$$
\begin{align}
g(v_0) = w_0 &\iff m_0 v_{0,x} + m_1 v_{0,y} + m_2 v_{0,z} + m_3 = w_0 \\
g(v_1) = w_1 &\iff m_0 v_{1,x} + m_1 v_{1,y} + m_2 v_{1,z} + m_3 = w_1 \\
g(v_2) = w_2 &\iff m_0 v_{2,x} + m_1 v_{2,y} + m_2 v_{2,z} + m_3 = w_2,
\end{align}
$$
which is a system of $3\times 3=9$ equations, but there are $3\times 4 = 12$ unknowns. This means that the system is underdetermined $\implies$ it has infinitely many solutions even if $v_0, v_1, v_2$ is a non-degenerate triangle (has non-zero area). There are various approaches to reduce the set of solutions to a single one.
Triangle Case: Solution
The first option is to supply a 3rd vector $e_3$ that is linearly independent of $e_1, e_2$ (doesn't lie in the same plane as $e_1$ and $e_2$), and similarly $f_3$ that isn't necessarily linearly independent of $f_1, f_2$. Then this becomes a special case of the tetrahedra transformations with $v_3 = v_0 + e_3$ and $w_3 = w_0 + f_3$. Depending on the direction and length of $e_3$ and $f_3$ the space outside of the planes in which the triangles lie gets sheared and scaled differently.
It also becomes geometrically clear at this point why we have $12$ unknowns and $9$ equations: the remaining $12-9=3$ unknowns correspond to this choice of 3-dimensional vectors. One would of course note that we had to specify two vectors: $e_3$ and $f_3$, which constitutes $6$ free variables instead of the $3$ pointed out above. The catch is that those $6$ variables get mixed together through $FE^{-1}$ in such a way that they could be represented through just $3$ free variables.
One simple option that guarantees linear independence is $e_3 = e_1 \times e_2$ and $f_3 = f_1 \times f_2$ where $\times$ is the cross product. The geometric implications of this choice are discussed in the TLDR: Differences Between Option 1 and Option 2 subsection. Another option would be to normalize the vectors: $e_3 = \frac{e_1 \times e_2}{\|e_1\times e_2\|}$ and $f_3 = \frac{f_1 \times f_2}{\|f_1\times f_2\|}$, where $\|\cdot\|$ is the Euclidean norm/length of the vectors. Finally, a choice of the form $e_3 = c (e_1\times e_2)$, where $c\ne0$, and $f_3=0$ results in a solution consistent with Option 2.
The difference to Option 1 is that points that are perpendicularly $X$ units away from the plane of the first triangle remain $X$ units away from the plane of the second triangle after transformation. That is, there is no dependence on the change of triangle area. Thus the case without normalization (i.e. Option 1) assumes the transformation is uniform rescaling with ratio $r = \frac{area(w_0,w_1,w_2)}{area(v_0,v_1,v_2)}$ -> rotation -> translation, while the normalized case assumes that scale may be applied only perpendicular to the normal of the triangle, followed by a rotation and translation.
Triangle Case: Pseudoinverse Solution
We take a slightly different approach by starting directly from the underdetermined system:
$$
\begin{align}
g(v_0) = w_0 &\iff m_0 v_{0,x} + m_1 v_{0,y} + m_2 v_{0,z} + m_3 = w_0 \\
g(v_1) = w_1 &\iff m_0 v_{1,x} + m_1 v_{1,y} + m_2 v_{1,z} + m_3 = w_1 \\
g(v_2) = w_2 &\iff m_0 v_{2,x} + m_1 v_{2,y} + m_2 v_{2,z} + m_3 = w_2.
\end{align}
$$
Rewriting $m_3 = w_0 - (m_0 v_{0,x} + m_1 v_{0,y} + m_2 v_{0,z})$ and plugging in the remaining two equations yields:
$$
\begin{gather}
m_0 e_{1,x} + m_1 e_{1,y} + m_2 e_{1,z} = f_1 \\
m_0 e_{2,x} + m_1 e_{2,y} + m_2 e_{2,z} = f_2 \\
\end{gather}
$$
We can write this more compactly in matrix form:
$$
\begin{gather}
M = \begin{bmatrix} m_0 & m_1 & m_2 \end{bmatrix}, \, E=\begin{bmatrix} e_1 & e_2 \end{bmatrix}, \, F=\begin{bmatrix} f_1 & f_2\end{bmatrix} \\
M E = F\\
ME = FI \\
ME = F(E^TE)^{-1}(E^TE) \\
ME = \left(F(E^TE)^{-1}E^T\right)E
\end{gather}
$$
In the above $(E^TE)^{-1}$ exists provided $e_1, e_2$ are linearly independent (the triangle $v_0,v_1,v_2$ is non-degenerate). By identifying the terms on the left- and right-hand sides we come to the conclusion that one possible form of $M$ is $M = F(E^TE)^{-1}E^T$.
The matrix $E^{+} = (E^TE)^{-1}E^T$ is in fact the Moore-Penrose pseudoinverse. It essentially projects points on the plane of the first triangle. For example, if we have a point $p$ then its projection is $p' = EE^+p$. Thus the final solution is:
$$
\begin{equation}
M = FE^{+}, \quad m_3 = w_0 - Mv_0.
\end{equation}
$$
Note the similarity to the solution $M = FE^{-1}$ where $E = \begin{bmatrix} e_1 & e_2 & e_3\end{bmatrix}$ and $F = \begin{bmatrix} f_1 & f_2 & f_3 \end{bmatrix}$. And actually the exact same solution for $M$ can be recovered by choosing $e_3 = c(e_1 \times e_2), c\ne 0$, and $f_3 = 0$.
For an efficient application of the Moore-Penrose pseudoinverse $E^{+}$ to some vector $p$ the following can be used:
$$
\begin{align}
E^{+}p &= (E^TE)^{-1}E^Tp \\
&= \left(\begin{bmatrix} e_1\cdot e_1 & e_1 \cdot e_2 \\ e_2 \cdot e_1 & e_2 \cdot e_2 \end{bmatrix}\right)^{-1}\begin{bmatrix} e_1^T \\ e_2^T \end{bmatrix}p \\
&= \frac{1}{(e_1\cdot e_1)(e_2 \cdot e_2) - (e_1\cdot e_2)^2}\begin{bmatrix} e_2\cdot e_2 & -e_1\cdot e_2 \\ -e_2 \cdot e_1 & e_1\cdot e_1 \end{bmatrix}\begin{bmatrix} e_1^T \\ e_2^T \end{bmatrix}p \\
&= \frac{1}{(e_1\cdot e_1)(e_2 \cdot e_2) - (e_1\cdot e_2)^2}\begin{bmatrix}(e_2\cdot e_2)(e_1\cdot p) - (e_1\cdot e_2)(e_2\cdot p) \\ (e_1\cdot e_1)(e_2\cdot p) -(e_1\cdot e_2)(e_1\cdot p)\end{bmatrix}.
\end{align}
$$
This results in 5 dot products, 6 scalar multiplications, and 3 subtractions for computing the application of $E^{+}$ to some vector $p$.
