derive (approximate) rotation transform matrix (numpy) on a unit sphere given a mapping of vectors (n=12)

While I am aware of vaguely similar questions, I seem to be stuck on this.

Buckminster Fuller introduced a spherical mapping of the world onto an icosahedron - it's known as the Dymaxion Map

A common way of identifying the cartesian coordinates of an icosahedron is by using the coordinates, where 𝜑 is the golden ratio: (1+√5)/2) or 2cos(π/5.0)

(0,±1,±1𝜑),(±1,±1𝜑,0),(±1𝜑,0,±1)

Expanding this out gives me the location of the 12 vertices of a regular icosahedron with side length 2:

    𝜑 = 2.0 * math.cos(π / 5.0)
    ico_vertices = [
        (0, -1, -𝜑), (0, -1, +𝜑), (0, +1, -𝜑), (0, +1, +𝜑),
        (-1, -𝜑, 0), (-1, +𝜑, 0), (+1, -𝜑, 0), (+1, +𝜑, 0),
        (-𝜑, 0, -1), (-𝜑, 0, +1), (+𝜑, 0, -1), (+𝜑, 0, +1)
    ]

Needless to say these need to be normalised.

    iv = np.array(ico_vertices)
    iv_n = ((iv[:, None] ** 2).sum(2) ** 0.5).reshape(-1, 1)
    ico = iv / iv_n  #This is the starting set of vertices.

Here is an image of the golden ratio 𝜑 coordinates projected onto Google Maps.

Fuller's icosahedron, mapped onto a spherical projection of the globe, defined using these 12 vertices: (the xyz values are already normalised to a unit sphere)

{
    "vertices":[
        {"name":"CHINA",                 "ll": [[39.10000000, "N"],[122.30000000,"E"]], "xyz":[-0.41468222, 0.65596241,  0.63067581]},
        {"name":"NORWAY",                "ll": [[64.70000000, "N"],[10.53619898 ,"E"]], "xyz":[ 0.42015243, 0.07814525,  0.90408255]},
        {"name":"ARABIAN SEA",           "ll": [[10.44734504, "N"],[58.15770555 ,"E"]], "xyz":[ 0.51883673, 0.83542038,  0.18133184]},
        {"name":"LIBERIA",               "ll": [[2.30088201 , "N"],[5.24539058  ,"W"]], "xyz":[ 0.99500944, -0.09134780, 0.04014717]},
        {"name":"PUERTO RICO",           "ll": [[23.71792533, "N"],[67.13232659 ,"W"]], "xyz":[ 0.35578140, -0.84358000, 0.40223423]},
        {"name":"ALASKA",                "ll": [[50.10320164, "N"],[143.47849033,"W"]], "xyz":[-0.51545596, -0.38171689, 0.76720099]},
        {"name":"BUENOS AIRES",          "ll": [[39.10000000, "S"],[57.70000000 ,"W"]], "xyz":[ 0.41468222, -0.65596241,-0.63067581]},
        {"name":"ANTARCTICA",            "ll": [[64.70000000, "S"],[169.46380102,"W"]], "xyz":[-0.42015243, -0.07814525,-0.90408255]},
        {"name":"PITCAIRN ISLAND",       "ll": [[10.44734504, "S"],[121.84229445,"W"]], "xyz":[-0.51883673, -0.83542038,-0.18133184]},
        {"name":"GILBERT ISLAND",        "ll": [[2.30088201 , "S"],[174.75460942,"E"]], "xyz":[-0.99500944, 0.09134780, -0.04014717]},
        {"name":"AUSTRALIA",             "ll": [[23.71792533, "S"],[112.86767341,"E"]], "xyz":[-0.35578140, 0.84358000, -0.40223423]},
        {"name":"PRINCE EDWARD ISLAND",  "ll": [[50.10320164, "S"],[36.52150967 ,"E"]], "xyz":[ 0.51545596, 0.38171689, -0.76720099]}
    ]
}

Here is an image of the dymaxion coordinates projected onto Google Maps.

The dymaxion coordinates are (via a json load) loaded into a numpy array 'dym_in'. The order of the two definitions is not the same - so the mapping is (this may be wrong).

    i2d = [6, 4, 10, 0, 8, 9, 3, 2, 7, 5, 11, 1]
    # ico[i] is equivalent to dym_in[i2d[i]]
    dym = np.array([dym_in[m] for m in i2d ])

So now I have 12 normalised vertices in 'ico' and 12 dymaxion map vertices in 'dym', which are ordered such that ico[x] => dym[x]. I want to find the rotation (or approximate rotation) matrix that transforms ico to dym.

I say approximate, because the coordinates in given dym may not exactly mathematically define an icosahedron. I do not know because I do not know how to derive the transform!

What I know for sure is that the geoid is not relevant here - the Dymaxion starts from a spherical earth projection.

Likewise, I freely admit there may be bugs in my assumptions above.

What I want is to be able to derive the rotational matrix of any set of 12 icosahedral points from the initial golden-ratio starting set - bearing in mind that there are several 12(?) rotations to choose from, of course.

Note: the permutation matrix i2d supplied in the original post is not correct and does not give a correct mapping of points (so NO method would be able to compute the rotation matrix from it). An alternative i2d array was found by searching permutations and is included in the code below. Note that, due to the symmetries of the icosahedron, there are many possible permutation arrays and this is just one. Note also that you do have to check that you have a pure rotation (det(R)=1) and not one including a reflection (det(R)=-1). Having just been caught out by that I've now put a determinant check in at the end.

You are asking for a 3x3 rotation matrix R taking position vectors U₁ to V₁, U₂ to V₂, U₃ to V₃ etc.

Choose any linearly independent triplet U₁, U₂, U₃. Then, as a matrix equation

(Make sure that the matrices are formed from successive column vectors.) Then just post-multiply by the inverse of the U matrix. This gives you R:

The code below produces the rotation matrix (for this particular mapping of icosahedral vertices). It also does a determinant check to make sure that you have a pure rotation and not an additional reflection (which could also map the point).

import math
import numpy as np

phi = 2.0 * math.cos(math.pi / 5.0)
ico = np.array( [ [   0,   -1, -phi],
                  [   0,   -1, +phi],
                  [   0,   +1, -phi],
                  [   0,   +1, +phi],
                  [  -1, -phi,    0],
                  [  -1, +phi,    0],
                  [  +1, -phi,    0],
                  [  +1, +phi,    0],
                  [-phi,    0,   -1],
                  [-phi,    0,   +1],
                  [+phi,    0,   -1],
                  [+phi,    0,   +1] ] )
for i in range( 12 ): ico[i] = ico[i] / np.linalg.norm( ico[i] )     # Normalise




dym_in = [
  [-0.41468222, 0.65596241,  0.63067581],
  [ 0.42015243, 0.07814525,  0.90408255],
  [ 0.51883673, 0.83542038,  0.18133184],
  [ 0.99500944, -0.09134780, 0.04014717],
  [ 0.35578140, -0.84358000, 0.40223423],
  [-0.51545596, -0.38171689, 0.76720099],
  [ 0.41468222, -0.65596241,-0.63067581],
  [-0.42015243, -0.07814525,-0.90408255],
  [-0.51883673, -0.83542038,-0.18133184],
  [-0.99500944, 0.09134780, -0.04014717],
  [-0.35578140, 0.84358000, -0.40223423],
  [ 0.51545596, 0.38171689, -0.76720099]
]

# ico[i] corresponds to dym_in[i2d[i]]
# i2d = [ 6, 4, 10, 0, 8, 9, 3, 2, 7, 5, 11, 1 ]    # Original permutation array: this is INCORRECT
i2d = [ 0, 3, 9, 6, 2, 7, 1, 8, 10, 11, 5, 4 ]      # Found by searching permutations
dym = np.array([dym_in[m] for m in i2d ])


U = np.zeros( ( 3, 3 ) )
V = np.zeros( ( 3, 3 ) )
independent = ( 0, 4, 8 )
for r in range( 3 ):
    for c in range( 3 ):
        U[r,c] = ico[independent[c],r]
        V[r,c] = dym[independent[c],r]

R = V @ np.linalg.inv( U )

for i in range( 12 ):
    print( "Vertex ", i, "   dym[i] = ", dym[i], "    R.ico = ", R @ (ico[i].T) )

print( "\nRotation matrix:\n", R )
print( "\nCheck determinant = ", np.linalg.det( R ) )

Output:

Vertex  0    dym[i] =  [-0.41468222  0.65596241  0.63067581]     R.ico =  [-0.41468222  0.65596241  0.63067581]
Vertex  1    dym[i] =  [ 0.99500944 -0.0913478   0.04014717]     R.ico =  [ 0.99500943 -0.0913478   0.04014718]
Vertex  2    dym[i] =  [-0.99500944  0.0913478  -0.04014717]     R.ico =  [-0.99500943  0.0913478  -0.04014718]
Vertex  3    dym[i] =  [ 0.41468222 -0.65596241 -0.63067581]     R.ico =  [ 0.41468222 -0.65596241 -0.63067581]
Vertex  4    dym[i] =  [0.51883673 0.83542038 0.18133184]     R.ico =  [0.51883673 0.83542038 0.18133184]
Vertex  5    dym[i] =  [-0.42015243 -0.07814525 -0.90408255]     R.ico =  [-0.42015243 -0.07814525 -0.90408256]
Vertex  6    dym[i] =  [0.42015243 0.07814525 0.90408255]     R.ico =  [0.42015243 0.07814525 0.90408256]
Vertex  7    dym[i] =  [-0.51883673 -0.83542038 -0.18133184]     R.ico =  [-0.51883673 -0.83542038 -0.18133184]
Vertex  8    dym[i] =  [-0.3557814   0.84358    -0.40223423]     R.ico =  [-0.3557814   0.84358    -0.40223423]
Vertex  9    dym[i] =  [ 0.51545596  0.38171689 -0.76720099]     R.ico =  [ 0.51545596  0.38171689 -0.767201  ]
Vertex  10    dym[i] =  [-0.51545596 -0.38171689  0.76720099]     R.ico =  [-0.51545596 -0.38171689  0.767201  ]
Vertex  11    dym[i] =  [ 0.3557814  -0.84358     0.40223423]     R.ico =  [ 0.3557814  -0.84358     0.40223423]

Rotation matrix:
 [[-0.09385435 -0.55192398  0.82859596]
 [-0.72021144 -0.53698041 -0.43925792]
 [ 0.68737678 -0.63799058 -0.34710402]]

Check determinant =  0.999999999723162