I have the following problem:
What are the vectors $\vec{r}_i$ that yield the maximal value of this determinant? \begin{equation} \det \left( \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \end{bmatrix} + \begin{bmatrix} \cdots & \vec{r}_1^T & \cdots \\ \cdots & \vec{r}_2^T & \cdots \\ \cdots & \vec{r}_3^T & \cdots \\ \cdots & \vec{r}_4^T & \cdots \end{bmatrix}_{4 \times 3} \begin{bmatrix} \vdots & \vdots & \vdots & \vdots \\ \vec{r}_1 & \vec{r}_2 & \vec{r}_3 & \vec{r}_4 \\ \vdots & \vdots & \vdots & \vdots \end{bmatrix}_{3 \times 4} \right) \end{equation} where $\vec{r}_i$ are all vectors in $\mathbb{R}^3$, s.t. $\left\vert \vec{r}_i \right\vert \leq 1 $. Simply put, the problem is maximizing the determinant of a matrix $A$ with $A_{ij} = 1+\vec{r}_i \cdot \vec{r}_j$, where $\vec{r}_i$ are all inside the unit sphere in three dimensions.
I think the solution is as follows: $\vec{r}_i$ all have unit norm, and if their end points form a tetrahedron. I.e.: \begin{equation} \vec{r}_1 = \begin{bmatrix} \sqrt{ 8/9 } \\ 0 \\ -1 / 3 \end{bmatrix} , \vec{r}_2 = \begin{bmatrix} -\sqrt{ 2/9 } \\ \sqrt{ 2/3 } \\ -1 / 3 \end{bmatrix} , \vec{r}_3 = \begin{bmatrix} -\sqrt{ 2/9 } \\ -\sqrt{ 2/3 } \\ -1 / 3 \end{bmatrix} , \vec{r}_4 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \end{equation} but I have no idea how to prove it. The only thing even remotely related I found when looking up, was the Cayley-Menger determinant: https://en.wikipedia.org/wiki/Cayley%E2%80%93Menger_determinant
and yet I couldn't find a clear connection between that and my problem.
Your determinant is $\det(X^TX)$, where $$ X=\pmatrix{1&1&1&1\\ \vec{r}_1&\vec{r}_2&\vec{r}_3&\vec{r}_4}. $$ Let $$ U=\pmatrix{\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\ \vec{r}_1&\vec{r}_2&\vec{r}_3&\vec{r}_4}. $$ Then $\det(X)=\sqrt{3}\det(U)$ and hence your determinant is equal to $3\det(U)^2$.
Since $\|\vec{r}_j\|\le1$, each column of $U$ has norm $\le\sqrt{\frac13+1}=\sqrt{\frac43}$. Therefore, by QR factorisation, $|\det(U)|\le\left(\sqrt{\frac43}\right)^4=\frac{16}{9}$ and equality holds if and only if $U$ is a scalar multiple of an orthogonal matrix. It follows that the maximum possible value of $3\det(U)^2$ is $3\left(\frac{16}{9}\right)^2=\frac{256}{27}$. This occurs, e.g. when $$ U=\pmatrix{\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\ \vec{r}_1&\vec{r}_2&\vec{r}_3&\vec{r}_4} =\frac{1}{\sqrt{3}}\pmatrix{1&1&1&1\\ 1&1&-1&-1\\ 1&-1&1&-1\\ 1&-1&-1&1}. $$