What you want is
$${\tiny{\begin{pmatrix}\textbf{v}_{1,x}\\\textbf{v}_{1,y}\\\textbf{v}_{1,z}\\\vdots\\\textbf{v}_{n,x}\\\textbf{v}_{n,y}\\\textbf{v}_{n,z}\end{pmatrix}}={\begin{pmatrix}\binom{N}{0}t_1^0(t_1^{N-0})& 0&0&\binom{N}{1}t_1^1(t_1^{N-1})&\ldots&0&\binom{N}{N}t_1^N(t_1^{N-N})&0&0
\\
0&\binom{N}{0}t_1^0(t_1^{N-0})& 0&0&\ldots&0&0&\binom{N}{N}t_1^N(t_1^{N-N})&0
\\
0&0&\binom{N}{0}t_1^0(t_1^{N-0})& 0&\ldots&\binom{N}{N-1}t_1^0(t_1^{N-N+1})&0&0&\binom{N}{N}t_1^N(t_1^{N-N})
\\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots
\\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots
\\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\
\binom{N}{0}t_n^0(t_n^{N-0})& 0&0&\binom{N}{1}t_n^1(t_n^{N-1})&\ldots&0&\binom{N}{N}t_n^N(t_n^{N-N})&0&0
\\
0&\binom{N}{0}t_n^0(t_n^{N-0})& 0&0&\ldots&0&0&\binom{N}{N}t_n^N(t_n^{N-N})&0
\\
0&0&\binom{N}{0}t_n^0(t_n^{N-0})& 0&\ldots&\binom{N}{N-1}t_n^0(t_n^{N-N+1})&0&0&\binom{N}{N}t_n^N(t_n^{N-N})\end{pmatrix}}{\begin{pmatrix}\textbf{P}_{0,x}\\\textbf{P}_{0,y}\\\textbf{P}_{0,z}\\\vdots\\\textbf{P}_{N,x}\\\textbf{P}_{N,y}\\\textbf{P}_{N,z}\end{pmatrix}}}$$
The matrix is the design matrix needed for the least square parameter:
$\textbf{P}_{vec}=(M^TM)^{-1}M^T\textbf{v}_{vec}$
See https://en.wikipedia.org/wiki/Ordinary_least_squares
Note when fitting a bezier, it usually doesn't matter if your specified points denoted by $\textbf{v}$ occurs exactly at those $\textbf{t}$ values you gave. The real fun begins when you let the $t$-values vary only preserving the order. If that is done you get a closer fit than the ordinary least squares.
Even more fun will arise if you want to minimize an other error than the least square error.