# Create a set of vectors from given pairwise dot product

Assume I am given a list of all the pairwise dot products of an unknown set of $$n$$ vectors. How can I efficiently generate such a set (or all possible sets) with the extra constraint that all vectors are unit vectors.

What I have tried

We can always choose the first vector to be the one that has 1 as its first component and 0 all the rest. Similarly, the second can have as the first element $$p_{12}$$ and the second $$\sqrt{1-p_{12}^2}$$, where $$p_{ij}$$ denotes the overlap between states $$i,j$$ and similarly for the rest. That is, $$v_1=e_1 \notag \\ v_2=p_{12}e_1 + \sqrt{1-p_{12}^2} e_2 \notag \\ \vdots$$ where $$e_i$$ denotes the $$i-$$th vector of an orthonormal basis.

My code for doing this is the following:

createStatesFromOverlaps[overlapMat_] :=
Module[{st, sol, dim, cf, e},
dim = Dimensions[overlapMat][[1]];
e[x_] := UnitVector[dim, x];
st = Table[
Sum[e[i]*If[i == 1, overlapMat[[i, x]], cf[x, i]], {i, 1, x}], {x,
1, dim}];
sol = Solve[
Table[
st[[i]] . st[[j]] == overlapMat[[i, j]], {i, 1, dim}, {j, 1,
i}] // Flatten,
Table[cf[i, j], {i, 1, dim}, {j, 2, i}] // Flatten];
Return[{st /. sol}]
]


The overlapMat is the given matrix $$P$$ of all possible overlaps, $$p_{ij}$$.

The issue

For two vectors this works well and it gives the correct result quickly usually. However, even for some simple 3 vector cases this takes ages to compute if the overlap matrix is exact and not numeric. If it is numeric, then it gives a warning about using Solve with inexact coefficients but otherwise gets a solution quickly. I was a bit surprised by that and was expecting the 3 vector case to be very quick, so I suspect that my implementation is not very clever.

Question: My question is, can I modify the code in some clever way to solve this problem efficiently?

Specific Example: Take as the matrix of overlaps the following one:

M={{1, ((1500 Sqrt[119] + 49 Sqrt[1365] + 7 Sqrt[23205] +
21 Sqrt[3885] (7 + Sqrt[17])) (3000 Sqrt[28886] +
Sqrt[30] (7 + Sqrt[17]) (111 + Sqrt[481])))/35000000000, (
Sqrt[3] (7 + Sqrt[17]) (Sqrt[65] + 3 Sqrt[185]))/
1000}, {((1500 Sqrt[119] + 49 Sqrt[1365] + 7 Sqrt[23205] +
21 Sqrt[3885] (7 + Sqrt[17])) (3000 Sqrt[28886] +
Sqrt[30] (7 + Sqrt[17]) (111 + Sqrt[481])))/35000000000, 1, (
3 (500 Sqrt[21] + 21 Sqrt[455] + 63 Sqrt[1295] + 3 Sqrt[7735] +
9 Sqrt[22015]) (200 Sqrt[
11110] + (7 + Sqrt[17]) (Sqrt[78] + 3 Sqrt[222])))/
1400000000}, {(Sqrt[3] (7 + Sqrt[17]) (Sqrt[65] + 3 Sqrt[185]))/
1000, (
3 (500 Sqrt[21] + 21 Sqrt[455] + 63 Sqrt[1295] + 3 Sqrt[7735] +
9 Sqrt[22015]) (200 Sqrt[
11110] + (7 + Sqrt[17]) (Sqrt[78] + 3 Sqrt[222])))/1400000000,
1}}


This seems to not compute in an expected time frame in my computer. However, it gives an immediate answer if I call the function as createStatesFromOverlaps[M//N] instead (but with the warning for inexact coefficients).

Thank you!

## 1 Answer

The solution is contained in $$U\sqrt V$$ where $$UVW^\star=X$$ is the singular value decomposition of your dot matrix $$X$$.

Remove["Global*"]
dim = 6;
numvecs = 20;
originalVectors = RandomPoint[Sphere[dim], numvecs];
dotmtx = dotm[originalVectors];

dotm[vs_?(MatrixQ[#, NumericQ] &)] := Outer[Dot, vs, vs, 1]
loss[vs_?(MatrixQ[#, NumericQ] &)] := Total[Abs[dotmtx - dotm[vs]], 2]

{u, v, w} = SingularValueDecomposition[dotmtx];
isovecs = u . MatrixPower[v, 1/2][[All, ;; dim]];


SingularValueDecomposition also works on symbolic matrices too. If we use your example:

Remove["Global*"]
dim = 3;
numvecs = 3;
originalVectors = RandomPoint[Sphere[dim], numvecs];
dotmtx = {{1, ((1500 Sqrt[119] + 49 Sqrt[1365] + 7 Sqrt[23205] +
21 Sqrt[3885] (7 + Sqrt[17])) (3000 Sqrt[28886] +
Sqrt[30] (7 + Sqrt[17]) (111 + Sqrt[481])))/
35000000000, (Sqrt[3] (7 + Sqrt[17]) (Sqrt[65] + 3 Sqrt[185]))/
1000}, {((1500 Sqrt[119] + 49 Sqrt[1365] + 7 Sqrt[23205] +
21 Sqrt[3885] (7 + Sqrt[17])) (3000 Sqrt[28886] +
Sqrt[30] (7 + Sqrt[17]) (111 + Sqrt[481])))/35000000000,
1, (3 (500 Sqrt[21] + 21 Sqrt[455] + 63 Sqrt[1295] +
3 Sqrt[7735] +
9 Sqrt[22015]) (200 Sqrt[
11110] + (7 + Sqrt[17]) (Sqrt[78] + 3 Sqrt[222])))/
1400000000}, {(Sqrt[3] (7 + Sqrt[17]) (Sqrt[65] + 3 Sqrt[185]))/
1000, (3 (500 Sqrt[21] + 21 Sqrt[455] + 63 Sqrt[1295] +
3 Sqrt[7735] +
9 Sqrt[22015]) (200 Sqrt[
11110] + (7 + Sqrt[17]) (Sqrt[78] + 3 Sqrt[222])))/
1400000000, 1}};

dotm[vs_?(MatrixQ[#, NumericQ] &)] := Outer[Dot, vs, vs, 1]
loss[vs_?(MatrixQ[#, NumericQ] &)] := Total[Abs[dotmtx - dotm[vs]], 2]

{u, v, w} = SingularValueDecomposition[dotmtx];
isovecs = u . MatrixPower[v, 1/2][[All, ;; dim]];
Norm /@ N[isovecs]
loss[N@isovecs]


The result isovecs is very complicated but numerically comes out as this (up to isometry):

{{0.978501, 0.154241, -0.136916}, {0.637032, -0.770185, 0.0317041}, {0.921297, 0.368728, 0.123495}}

• That is really nice! It can calculate the vectors exactly quickly (even if the answer is pretty horrible). That is what I was looking for. Thank you Jan 28, 2023 at 3:00