title: efficiently calculate list of 3d rotation matrices in numpy or scipy
body: <p>I have a list of N unit-normalized 3D vectors <strong>p</strong> stored in a numpy ndarray with shape (N, 3). I have another such list, <strong>q</strong>. I want to calculate an ndarray <strong>U</strong> of shape (N, 3, 3) storing the rotation matrices that rotate each point in <strong>p</strong> to the corresponding point <strong>q</strong>. </p>
  
  <p>The list of rotation matrices U should satisfy:</p>
  
  <pre><code>np.all(np.einsum('ijk,ik-&gt;ij', U, p) ==  q)
  </code></pre>
  
  <p>On a point-by-point basis, the problem reduces to being able to compute a rotation matrix for a rotation of some angle about some axis. Code solving the single-point case appears below:</p>
  
  <pre><code>def rotation_matrix(angle, direction):
      direction = np.atleast_1d(direction).astype('f4')
      sina = np.sin(angle)
      cosa = np.cos(angle)
      direction = direction/np.sqrt(np.sum(direction*direction))
  
      R = np.diag([cosa, cosa, cosa])
      R += np.outer(direction, direction) * (1.0 - cosa)
      direction *= sina
      R += np.array(((0.0, -direction[2], direction[1]),
                  (direction[2], 0.0, -direction[0]),
                  (-direction[1], direction[0],  0.0)))
      return R
  </code></pre>
  
  <p>What I need is a function that behaves exactly as the above function, but instead of accepting a single angle and a single direction, it accepts an <code>angles</code> array of shape (npts, ) and a <code>directions</code> array of shape (npts, 3). The code below is only partially finished - the problem is that neither <code>np.diag</code> nor <code>np.outer</code> accept an <code>axis</code> argument </p>
  
  <pre><code>def rotation_matrices(angles, directions):
      directions = np.atleast_2d(directions)
      angles = np.atleast_1d(angles)
      npts = directions.shape[0]
      directions = directions/np.sqrt(np.sum(directions*directions, axis=1)).reshape((npts, 1))
  
      sina = np.sin(angles)
      cosa = np.cos(angles)
  
      #  Lines below require extension to 2d case - np.diag and np.outer do not support axis arguments
      R = np.diag([cosa, cosa, cosa])
      R += np.outer(directions, directions) * (1.0 - cosa)
      directions *= sina
      R += np.array(((0.0, -directions[2], directions[1]),
                  (directions[2], 0.0, -directions[0]),
                  (-directions[1], directions[0],  0.0)))
      return R
  </code></pre>
  
  <p>Does either numpy or scipy have a compact vectorized function computing the appropriate rotation matrices in a way that avoids using for loops? The problem is that neither <code>np.diag</code> nor <code>np.outer</code> accept <code>axis</code> as an argument. My application will have N be very large, 1e7 or greater, so a vectorized function that keeps all the relevant axes aligned is necessary for performance reasons. </p>