# How to linearize Quaternions?

Based on an answer to one of my questions and the comments exchanged here earlier I noticed that I cannot uniformly sample Quaternion vectors for rotation even though if I know the bounds of each element of the vector. More specifically, I mean this comment:

If you're going to do your max/min approach to find bounds, I'd recommend that you: 1) throw out quaternions with small L2 norms (say, less than 12) as well as large ones, as otherwise your bounds won't be useful, and 2) generate a really huge number of quaternions to find your bounds (leave your computer on overnight), as otherwise you might miss some low-probability part of the region

To me, this means that Quaternions vectors of rotation are not linear and there they have some sort of weird curvature around the corners (or maybe in the middle too?).

So I wonder, is it possible to develop a method to linearize Quaternions?

The reasons that I want to linearize Quaternion vectors of rotation is as follow:

• I want to $$\pm$$ some value to each element of the Quaternion vector to avoid some degeneration issues that the optimization would face. If I simply add some values to a non-linearized Quaternion vector then I might drastically move in the space, which is not what I want (this is the main reason)
• I'm trying to sample Quaternion rotation vectors uniformly and without doing rejection sampling, which is computationally intensive on large-scale.
• I am using an optimization procedure to find some Quaternion rotation vector that fits my data. Sampling from Euler angles and then converting it to Quaternion makes the optimization procedure way harder because Euler angles are non-linear and hard to interpret.

Ideally, I would like to get a deterministic, one-to-one mapping function for each element of the Quaternion vector $$(x, y, z, w)$$ but I am also open to use some computational approaches to empirically learn/compute this linear mapping.

• What is your primary goal here? If you're trying to uniformly and randomly sample rotations, there's a straightforward way to do that that doesn't involve 'linearization' in any meaningful way... – Steven Stadnicki Nov 9 '18 at 22:29
• @StevenStadnicki I'm trying to sample Quaternion rotation vectors uniformly and without doing rejection sampling. I thought this is not possible based on the answer and comments to my earlier question. – Amir Nov 9 '18 at 22:32
• @JunekeyJeon That doesn't actually correctly sample uniformly; the distribution over the rotation angle isn't uniform the way it is in the case one dimension down. See, e.g., en.wikipedia.org/wiki/… – Steven Stadnicki Nov 9 '18 at 22:37
• You're going to have to compromise somewhere, because rotations are inherently nonlinear - as exhibited, for instance, by the fact that they're noncommutative. Being able to find the most appropriate set of linear approximations for a particular application is a very hard problem and requires a lot of domain-specific knowledge; I'm pretty sure people have gotten PhD theses out of comparable problems. – Steven Stadnicki Nov 9 '18 at 23:05
• I often do exactly this to generate an uniformly random orientation, using the rejection method. It has about 50% reject ratio for this, and easily verified uniformity. A 50% reject ratio is not computationally intensive, or you are using the wrong pseudorandom number generator. – Nominal Animal Nov 10 '18 at 7:08

I think it's better to provide this as a separate answer, because the contents are rather largely independent to the previous answer.

As I understood, what you want to do is a kind of gradient descent inside the Lie group $$\mathrm{SO}(3)$$ (the set of 3-dimensional rotation matrices; recall that a $$3\times3$$-matrix $$R$$ is a rotation matrix if and only if $$RR^{T}=I$$ and $$\mathrm{det}(R)=1$$) or $$\mathrm{SU}(2)$$ (the set of unit quaternions). I've done before some sort of a Gauss-Newton optimization inside a Lie group, so I guess I can help you.

I am not sure how much you know about differential geometry, but I guess you are not very familiar with concepts such as Lie group or tangent space. Thus, I tried to avoid sophisticated mathematical constructions as much as possible.

First of all, I want to emphasize that as long as what you are interested in is the rotation of vectors in $$\mathbb{R}^{3}$$, then the exact way of representing rotation should not matter. That means, there should be no difference (aside from numerical issues like floating-point errors and computational performance) no matter how you represent rotations, because there are ways to convert between those representations, which implies that an algorithm for some representation should have a translation in terms of another representation.

Given this, it is helpful to conceptually think what you are finding is a rotation matrix $$R\in\mathrm{SO}(3)$$ (because that is "the most natural" representation of rotation), although you might use quaternions or Euler angles for actual computing.

Now, let $$f:\mathrm{SO}(3)\rightarrow[0,\infty)$$ be the cost function we wish to minimize using a gradient descent-like method. Let $$R_{0}\in\mathrm{SO}(3)$$ be the current (or initial) estimation of the optimal point. Of course, you can write $$f$$ as a function of $$9$$ entries of the rotation matrix $$R$$ and compute the gradient in terms of those components, but this is certainly a mistake, because then the updated matrix might not be a rotation matrix anymore. This is a problem of exactly the same kind of yours, right? So, what we need here is a way to represent $$\mathrm{SO}(3)$$ using independent parameters. (I guess this is what you wanted to do with a "linearization.") Mathematicians call such a parametrization a coordinate map, a coordinate chart, or a chart, depending on the literature.

There is a canonical way of doing this, indeed. For a given vector $$\omega\in\mathbb{R}^{3}$$, think of the anti-symmetric matrix $$[\omega]_{\times}:=\begin{bmatrix}0&-\omega_{z}&\omega_{y}\\\ \omega_{z}&0&-\omega_{x}\\\ -\omega_{y}&\omega_{x}&0\end{bmatrix}.$$ One can show that, the matrix exponential $$\exp([\omega]_{\times})=I+[\omega]_{\times}+\frac{[\omega]_{\times}^{2}}{2}+\ \cdots\ \stackrel{[1]}{=}I+\frac{\sin\|\omega\|}{\|\omega\|}[\omega]_{\times}+\frac{1-\cos\|\omega\|}{\|\omega\|^{2}}[\omega]_{\times}^{2}$$ (to see why [1] holds, see https://en.wikipedia.org/wiki/Axis%E2%80%93angle_representation, the section explaining the Exponential map)

is the rotation matrix along the axis $$\frac{\omega}{\|\omega\|}$$ with the angle $$\|\omega\|$$. (The $$\omega v$$ in the previous answer exactly is the $$\omega$$ here.) And conversely, any rotation can be written in this way. The three parameters $$(\omega_{x},\omega_{y},\omega_{z})$$ are completely independent, so can be freely set. It does not have any constraint to satisfy. Thus, taking the gradient of $$f$$ with respect to these variables and then making an update according to it, is a completely legal algorithm.

However, often this is not a good method. I can't explain the reason very clearly, but I think one can say that that's because the coordinate map $$(\omega_{x},\omega_{y},\omega_{z})\mapsto\exp([\omega]_{\times})$$ is "good" only near the identity. I think you can verify yourself that, usually the things become very complicated if you attempt to write the cost function $$f$$ in terms of $$(\omega_{x},\omega_{y},\omega_{z})$$. (Do you see those horrible functions $$\frac{\sin\|\omega\|}{\|\omega\|}$$ and $$\frac{1-\cos\|\omega\|}{\|\omega\|^{2}}$$?) The resulting gradient descent method usually performs really poorly; even naive renormalization method in the quaternion space often performs better. (By renormalization method I mean, (1) just compute the derivative of $$f$$ with respect to four components of the quaternion, (2) update the quaternion accordingly, and then (3) normalize the resulting quaternion. Repeating (1)-(3), although heuristic at best, is often a valid method computing the optimal point of $$f$$. I think this renormalization method is close to what you were thinking about.)

So, what should we do then? I said the coordinate map we have talked about is "good" only near the identity. Indeed, the coordinate map is really good enough at least near the identity. Therefore, all you need is just to "translate" the coordinate map to the point $$R_{0}$$ by multiplying $$R_{0}$$; that is, instead of $$(\omega_{x},\omega_{y},\omega_{z})\mapsto\exp([\omega]_{\times})$$, you use $$(\omega_{x},\omega_{y},\omega_{z})\mapsto\exp([\omega]_{\times})R_{0}$$. (Depending on the situation, you might need to use $$R_{0}\exp([\omega]_{\times})$$ instead of $$\exp([\omega]_{\times})R_{0}$$.)

Let me give you an example that can help understanding what's going on. Let's say $$f$$ is given as: $$f(R):=\frac{1}{2}\|q-Rp\|^{2}.$$ This means, the cost is the squared-distance between two prescribed points $$p$$ and $$q$$, after rotating the point $$p$$ by the given rotation matrix $$R$$. Thus, the optimal rotation should align $$q$$ and $$Rp$$ in the same direction. Well, such a rotation can be computed directly, but let us just use this simple example to illustrate how gradient descent method can be adopted to $$\mathrm{SO}(3)$$.

Note that the problem of minimizing $$f(R)$$ is equivalent to that of $$\tilde{f}(\omega):=f(\exp([\omega]_{\times})R_{0})$$. The latter is a function on $$\mathbb{R}^{3}$$, thus we can make an attempt to do gradient descent update. Note that with respect to $$\tilde{f}$$, our current estimate $$R_{0}$$ corresponds to the origin $$\omega=0$$.

Now, let us compute the gradient of $$\tilde{f}$$ at the origin. Actually, it is easier (and more natural) to compute the directional derivative of $$\tilde{f}$$ rather than the gradient. Recall from vector calculus that $$\partial_{\omega}\tilde{f}=\langle\nabla\tilde{f},\omega\rangle,$$ where by $$\partial_{\omega}\tilde{f}$$ we denote the directional derivative of $$\tilde{f}$$ along the direction $$\omega$$. Thus, we can get $$\nabla\tilde{f}$$ by finding the directional derivative $$\partial_{\omega}\tilde{f}$$ of $$\tilde{f}$$ along an arbitrary direction $$\omega$$. In particular, at the origin, $$\partial_{\omega}\tilde{f}(0)=\lim_{\epsilon\rightarrow0} \frac{\tilde{f}(\epsilon\omega)-\tilde{f}(0)}{\epsilon} =\lim_{\epsilon\rightarrow0}\frac{f(\exp(\epsilon[\omega]_{\times})R_{0})-f(R_{0})}{\epsilon}.$$

Note that \begin{align*} f(R)&=\frac{1}{2}\langle q-Rp,q-Rp\rangle\\ &=\frac{1}{2}\left(\|q\|^{2}-\langle Rp,q\rangle-\langle q,Rp\rangle+\|Rp\|^{2}\right)\\ &=\frac{\|p\|^{2}+\|q\|^{2}}{2}-\langle q,Rp\rangle. \end{align*} (Note that $$\|Rp\|=p$$; the length of a vector is invariant under rotation.)

Therefore, for a small $$\epsilon>0$$, we get $$\frac{f(\exp(\epsilon[\omega]_{\times})R_{0})-f(R_{0})}{\epsilon} =-\left\langle q,\frac{\exp(\epsilon[\omega]_{\times})-I}{\epsilon}R_{0}p\right\rangle.$$ We can see from the Taylor expansion of $$\exp(\epsilon[\omega]_{\times})$$ that the right-hand side goes to $$-\langle q,[\omega]_{\times}R_{0}p\rangle$$ as $$\epsilon\rightarrow0$$. It can be easily seen that for any $$v\in\mathbb{R}^{3}$$, $$[\omega]_{\times}v=\omega\times v$$, the cross product of $$\omega$$ and $$v$$. Hence, $$\partial_{\omega}\tilde{f}(0)=-\langle q,\omega\times R_{0}p\rangle.$$ Using the identity $$\langle a,b\times c\rangle=\langle b,c\times a\rangle$$, we can rewrite this as $$-\langle \omega,R_{0}p\times q\rangle.$$

Comparing this with $$\partial_{\omega}\tilde{f}=\langle\nabla\tilde{f},\omega\rangle,$$ we get the conclusion: $$\nabla\tilde{f}(0)=-R_{0}p\times q.$$

Let's say the rate of movement is $$\alpha>0$$. Then, the next estimate of the optimal point of $$\tilde{f}$$ is $$\alpha R_{0}p\times q$$. In other words, the next estimate of the optimal rotation is: $$R_{1}=\exp(\alpha[R_{0}p\times q]_{\times})R_{0}.$$ This is the gradient descent update formula of $$f$$. Then, we can iterate this process to get $$R_{2}=\exp(\alpha[R_{1}p\times q]_{\times})R_{1}$$, $$R_{3}=\exp(\alpha[R_{2}p\times q]_{\times})R_{2}$$, and so on. Note that at each iteration we effectively use a "different cost function," because at the $$k$$th iteration what we manipulate is the function $$\tilde{f}(\omega)=f(\exp([\omega]_{\times})R_{k-1})$$.

If we use quaternions to represent rotations, the formula $$R_{1}=\exp(\alpha[R_{0}p\times q]_{\times})R_{0}$$ becomes even simpler to be evaluated, because the unit quaternion corresponding to the rotation matrix $$\exp([\omega]_{\times})$$ is just $$q=\left(\frac{\omega}{\|\omega\|}\sin\frac{\|\omega\|}{2}, \cos\frac{\|\omega\|}{2}\right)$$ in the $$q=(x,y,z,w)$$ convention. (Please refer to the previous answer.)

It is quite common that the cost function $$f$$ is of the form $$f(R)=\sum_{i=1}^{N}f_{i}(Rp_{i})$$, where $$p_{i}\in\mathbb{R}^{3}$$ and $$f_{i}:\mathbb{R}^{3}\rightarrow\mathbb{R}$$. After some computation as we did above, one can show by using chain rule that the directional derivative of $$\tilde{f}(\omega):=f(\exp([\omega]_{\times})R_{0})$$ along $$\omega$$ evaluated at the origin is $$\sum_{i=1}^{N}\langle\nabla f_{i}(R_{0}p_{i}),\omega\times R_{0}p_{i}\rangle =\sum_{i=1}^{N}\langle\omega,R_{0}p_{i}\times\nabla f_{i}(R_{0}p_{i})\rangle.$$ Hence, $$\nabla\tilde{f}(0)=\sum_{i=1}^{N}\big(R_{0}p_{i}\times\nabla f_{i}(R_{0}p_{i})\big),$$ thus the update formula is $$R_{1}=\exp\left(-\alpha\left[\sum_{i=1}^{N}\big(R_{0}p_{i}\times\nabla f_{i}(R_{0}p_{i})\big)\right]_{\times}\right)R_{0}.$$

Thank you for reading this loooong answer. I hope this answer helped you, but I think it was not clear enough to be easily grasped. Or, it could be possible that what I explained has nothing to do with your actual problem. So, please feel free to ask further.

• I really appreciate you spending so much time trying to explain things to me in detail. I am not entirely familiar with some of the stuff mentioned but it looks like this should put me into the right direction and this is possible the answer to my question. – Amir Nov 13 '18 at 15:35
• @Amir You're welcome. I edited the answer to make it more comprehensible, mainly the example section. As I wrote in the answer, please feel free to ask further. – Junekey Jeon Nov 13 '18 at 16:39

You can generate a uniformly sampled rotation by sampling the axis and the angle separately. As pointed out by Steven, the angle should not be uniformly sampled, though. I guess you already know how to obtain a quaternion from the axis-angle representation.

The complete procedure is given as follows:

1. Rotation angle: uniformly sample in the interval $$[0,\pi]$$, and find $$\omega$$ such that $$\omega-\sin\omega$$ equals to that sampled value.

2. Longitude of the rotation axis: uniformly sample in the interval $$[0,2]$$, and find $$\theta$$ such that $$1-\cos\theta$$ equals to the sampled value.

3. Latitude of the rotation axis: this is easy; just uniformly sample $$\phi$$ in the interval $$[0,2\pi]$$.

4. Compute the quaternion $$q=(\cos\frac{\omega}{2},v\sin\frac{\omega}{2})$$, where $$v=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$$.

Step 1 might be nontrivial. The inverse function of $$\omega-\sin\omega$$ is not an elementary function, but is a well-known function. See https://en.wikipedia.org/wiki/Kepler%27s_equation, for example. I think Newton's method will converge very fast.

In Step 2, you actually don't need to find $$\theta$$, because all you need are $$\sin\theta$$ and $$\cos\theta$$. Since $$\sin\theta$$ is guaranteed to be nonnegative, you can apply the formula $$\sin\theta=\sqrt{1-\cos^{2}\theta}$$.

There is another method which is not very pretty but usually practically better than the method above: sample four values randomly and independently according to the standard Gausssian distribution. Since the isotropic $$n$$-dimensional Gaussian distribution is just the $$n$$-fold product of 1-dimensional Gaussian distributions, the distribution of the obtained 4d vector has rotational symmetry. Thus you just need to normalize the resulting vector, then you get a uniformly sampled unit quaternion. Gaussian sampling is a well-known problem, so there are many known solutions. But to be precise, you might need to handle the case when the obtained vector is too small or too large. Then it becomes basically a sort of rejection sampling, which you wanted to avoid.

EDIT

I just found out in the Internet that in fact sampling Euler angles is a far easier way of getting unbiased, uniformly random rotation.

There are hundreds of different conventions of "Euler angles." The one used here is the following: $$R(\theta,\varphi_{1},\varphi_{2}):=Z(\varphi_{2})X(\theta)Z(\varphi_{1})$$ where $$Z(\varphi)=\begin{bmatrix}\cos\varphi&-\sin\varphi&0\\\sin\varphi&\cos\varphi&0\\0&0&1\end{bmatrix},\ X(\theta)=\begin{bmatrix}1&0&0\\0&\cos\theta&-\sin\theta\\0&\sin\theta&\cos\theta\end{bmatrix}.$$

The complete procedure is given as following:

1. Randomly sample $$\varphi_{1},\varphi_{2}$$ in $$[0,2\pi]$$.
2. Randomly sample $$u\in[0,2]$$, and find $$\theta\in[0,\pi]$$ such that $$1-\cos\theta=u$$; thus, $$\theta=\cos^{-1}(1-u)$$. However, in fact you don't need to compute arccos, becauses all you need are $$\sin\theta$$ and $$\cos\theta$$. In the range $$[0,\pi]$$, $$\sin\theta$$ is nonnegative, so $$\cos\theta=1-u$$ and $$\sin\theta=\sqrt{1-(1-u)^{2}}$$.
3. Compute the rotation matrix $$R(\theta,\varphi_{1},\varphi_{2})$$ as the above.
4. If you need a quaternion from the rotation matrix, then use the formula given in, for example, https://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Rotation_matrix_%E2%86%94_quaternion.

I can't easily explain why this gives a uniformly randomly sampled rotation. In a nutshell, that's because "$$\sin\theta\,d\theta\,d\varphi_{1}\,d\varphi_{2}$$ is a Haar measure of $$\mathrm{SO}(3)$$". I think you can safely consider "Haar measure" just as a different name of uniform distribution. But how to verify? It's a bit involved, and honestly I didn't check it myself. I referred to an online note http://www.math.ubc.ca/~feldman/m606/haar.pdf. And this paper https://www.hindawi.com/journals/tsm/1988/595067/abs/ also states the same claim :the equation (2). (You can regard what I used before the edit is the equation (4))

Q1) I don't know what you mean by sampling the axis and angle separately. Do you mean sampling Euler angle rotation vectors (rx, ry, rz)?

• No, they are completely different things. (And I guess the convention of Euler angles you are using is also very different from that used here.) Axis-angle representation is arguably the "most intuitive" way of representing a rotation. You specify an axis ($$v$$, a unit vector in 3D space), and specify how much angle $$\omega$$ you want to rotate along that axis.

Q2) Could you please give a more detailed example? I am a bit confused on how you construct 4 numbers at the end. At step 4, you say "Compute the quaternion" but I don't see how you get the 4 numbers at the end. Quaternion rotation vectors are defined as (rx, ry, rz, rw) for me. Are you thinking of a different representation?

• The quaternion corresponding to $$\omega$$ and $$v$$ can be obtained as $$q=(v_{x}\sin\frac{\omega}{2},v_{y}\sin\frac{\omega}{2},v_{z}\sin\frac{\omega}{2},\cos\frac{\omega}{2})$$, where $$v=(v_{x},v_{y},v_{z})$$ is the (normalized) axis of rotation. Here, I used the convention $$q=(x,y,z,w)$$ because that seems to be the one you are using; before the edit, I used the convention $$q=(w,x,y,z)$$.

Q3) In case of sampling from a Normal distribution, how should I measure how small or large the vector is? Should I calculate the L2 distance? If so, what would be good numbers for determining whether the L2 distance is too large or small so that I reject that sample?

• Yes, you should use the $$L^{2}$$-norm to normalize the 4D vector. There is no simple and perfect way of determining the threshold of rejection, because that largely depend on the situation. For example, the "optimal answer" might depend on whether you use float or double, whether you use precise/fast floating-point computation modes, the scaling of the Gaussian distribution you are using, etc.. Just you need to ensure that $$r^{2}=x^{2}+y^{2}+z^{2}+w^{2}$$ is an "ordinary" floating point number so that you don't run into a trouble computing $$(\frac{x}{r},\frac{y}{r},\frac{z}{r},\frac{w}{r})$$. There are tricks that help to mitigate this issue. For example, pre-dividing $$(x,y,z,w)$$ by $$\max\{|x|,|y|,|z|,|w|\}$$ is a good way to resolve lots of numerical troubles. Of course, however, doing so makes your computation slower. Anyway, I think the new method seems to be easy enough to be better than this Gaussian method.
• Thank you for the answer and explanations. Just one thing ... . I am not good at rotation and Quaternion as much as you probably imagined I am. So I want to ask you a couple of question to make things more clear. If possible, I would appreciate if you can edit you answer if you feel you need to so that I can understand things easier. Q1) I don't know what you mean by sampling the axis and angle separately. Do you mean sampling Euler angle rotation vectors (rx, ry, rz)? – Amir Nov 10 '18 at 3:26
• Q2) Could you please give a more detailed example? I am a bit confused on how you construct 4 numbers at the end. At step 4, you say "Compute the quaternion" but I don't see how you get the 4 numbers at the end. Quaternion rotation vectors are defined as (rx, ry, rz, rw) for me. Are you thinking of a different representation? Q3) In case of sampling from a Normal distribution, how should I measure how small or large the vector is?Should I calculate the $L_2$ distance? If so,what would be good numbers for determining whether the $L_2$ distance is too large or small so that I reject that sample? – Amir Nov 10 '18 at 3:31
• Thanks, things are more clear. However,I still have not got an answer to my question.My point is to be able to interpolate between elements of the Quaternion vector meaning that if the bound is $[0, 1]$ for each of the elements and I add $0.1$ to all four elements, the rotation is only changed by $10%$. I couldn't come up with a better explanation of what I want and maybe saying something like "changing the rotation by $10%$" would be a bit vague but I hope I am able to better convey what I want. Do you have any ideas on how to do this? Is there any computational way to linearize Quaternions? – Amir Nov 10 '18 at 17:09
• @Amir Before diving into detailed research, I think I need to understand your intent clearer. Q1) First, I thought you were doing some kind of Monte Carlo method and that's the reason why you need a random quaternion. Now, I think you aren't, and perhaps you were solving a kind of optimization problem. Is it right? Q2) I guess the reason why you want to "perturb" a given quaternion is because you are doing some kind of gradient descent or related optimization procedure. And I guess you are worrying if the updated quaternion might become something invalid. Am I correct? – Junekey Jeon Nov 11 '18 at 12:20
• Thanks for trying to clarify. The answer is yes to both of your questions. I was partly interested in the sampling Quaternions partly because I was thinking it can potentially resolves the possibility of getting "invalid" gradient descent updates and eventually help the optimization. – Amir Nov 11 '18 at 15:12

IMO a suitable solution is (as others said) by working with the axis/amplitude rotation system (Rodrigues form). Sample the amplitude uniformly. Sampling the axis direction uniformly is more tricky. You can work this out in spherical coordinates. Sample the azimuth uniformly, and the elevation non-uniformly, noting that the length of a parallel is proportional to the sine of the elevation. (If I am right, taking the arc cosine of a uniform variable does the trick.)

Also check https://en.wikipedia.org/wiki/Lambert_cylindrical_equal-area_projection for a uniform point density on a sphere.

Notice that it is impossible to avoid singularities.

• Thank you but the main problem is not sampling. I still have not got an answer to my question .My point is to be able to interpolate between elements of the Quaternion vector meaning that if the bound is $[0, 1]$ for each of the elements and I add $0.1$ to all four elements, the rotation is only changed by $10%$.I couldn't come up with a better explanation of what I want and maybe saying something like "changing the rotation by $10%$" would be vague but I hope I am able to better convey what I want.Do you have any ideas on how to do this?Is there any computational way to linearize Quaternions? – Amir Nov 10 '18 at 17:10
• @Amir: My answer implicitly says it. You can interpolate linearly on the amplitude, the azimuth and the argument of the arc cosine. – Yves Daoust Nov 11 '18 at 9:45

In a comment, OP was surprised that a rejection method implementation with a rejection ratio of about 50% is possible. So, I'll outline the rejection method based approach here, even though OP ruled that out in the question. Therefore, consider this an extended comment and not an answer proper.

Unit quaternions or versors describing an orientation are tightly coupled to the axis-angle representation. If $$\hat{a}$$ is an unit vector that describes the direction of the rotation axis, and $$\theta$$ is the rotation around that axis, then versor (unit quaternion) $$\mathbf{q}$$ is $$\bbox{\mathbf{q} = \cos\left(\frac{\theta}{2}\right) + \sin\left(\frac{\theta}{2}\right) \hat{a}}$$

To get an uniform distribution of orientations, we only need to generate an axis vector in an uniform direction, and an angle with an uniform distribution over any range of angles a positive integer multiple of $$180\text{°}$$.

Let's look at how this is implemented in pseudocode. Let's assume you have a function Prnd() that yields an uniform random number in $$[0, 1)$$, and that Sin() and Cos() take the angle in radians:

Function random_versor():

Do:
Let  x = 2*Prnd() - 1    # [-1, +1)
Let  y = 2*Prnd() - 1    # [-1, +1)
Let  z = 2*Prnd() - 1    # [-1, +1)
Let  nn = x*x + y*y + z*z
While (nn < 0.25) Or (nn >= 1.0)

Let theta = Prnd() * π * 0.5

Let scale = Sin(theta) / n

q.w = Cos(theta)
q.x = scale * x
q.y = scale * y
q.z = scale * z

Return q
End Function


We generate the initial vectors in the cube $$(-1..+1, -1..+1, -1..+1)$$, which has volume $$2^3 = 8$$. We only accept the ones that are within the unit sphere, but further than $$\sqrt{0.25} = 0.5$$ from the center; volume $$4π/3 - 4π(0.5)^3/3 \approx 3.665$$. The ratio is $$3.665/8 \approx 46 \%$$, so on average, about half the candidate random points are rejected.

Because the volume is spherical, treating each point as a vector, and scaling them to unit length, produces an uniform distribution of points on the shell of the unit sphere; and therefore uniform random unit vectors.

There is a method of generating such unit vectors directly, of course: $$\bbox{\begin{cases} \theta = 2 \pi u , & 0 \le u \le 1 \\ \phi = \arccos(2 v - 1) , & 0 \le v \le 1 \\ \end{cases}}$$ where $$u$$ and $$v$$ are uniform random numbers between zero and one; and $$\bbox{\begin{cases} x = \sin(\phi) \cos(\theta) \\ y = \sin(\phi) \sin(\theta) \\ z = \cos(\phi) \\ \end{cases}}$$

In pseudocode, that corresponds to

Function random_versor():

Let  u = Prnd()   # [0, 1)
Let  v = Prnd()   # [0, 1)

Let  theta = 2 * π * u
Let  phi = ArcCos(2*v - 1)
Let  beta = Prnd() * π * 0.5   # Half the angle around the axis

Let  b = Sin(beta)
Let  s = b * Sin(phi)

q.w = Cos(beta)
w.x = s * Cos(theta)
w.y = s * Sin(theta)
w.z = b * Cos(phi)

Return q
End Function


In practice, the cost of the average extra three Prng() calls and the one square root call per versor generated is less than the cost of the $$\arccos$$ and extra $$\sin$$ and $$\cos$$, which means that for sane pseudorandom number generators like Mersenne Twister and Xorshift variants, the exclusion method is more efficient.

To verify either generator, use 0 for the rotation around the axis, and check if the generated axis vectors are uniformly distributed on the unit sphere.

Because of the relationship between the $$w$$ and the $$x$$, $$y$$, $$z$$ components, it is not possible for them to be generated independently, and still have the result be normalizable to an unit quaternion describing orientations with an uniform distribution. The "linearization" OP describes, is impossible.

Do remember that for $$\mathbf{q} = (w, x, y, z)$$ to be an unit quaternion AKA versor, $$w^2 + x^2 + y^2 + z^2 = 1$$ but that it isn't just a four-component vector; it describes a rotation of $$2\arccos(w)$$ around the axis in direction $$(x, y, z)$$.

In the other question, it seems to me OP is trying to find a way to "perturb" an orientation slightly, in a random fashion.

In my opinion, the easiest way to do that, is to use the unit axis $$\hat{a}$$ and rotation $$\theta$$ used to create the versor, and perturb those, rather than trying to perturb the versor directly.

You can interpolate between any two versors $$\mathbf{q}_0$$ and $$\mathbf{q}_1$$ using $$\bbox{\mathbf{q} = \frac{ (1 - \lambda) \mathbf{q}_0 + \lambda \mathbf{q}_1 }{\left\lVert (1 - \lambda) \mathbf{q}_0 + \lambda \mathbf{q}_1 \right\rVert}}$$ where as usual, $$\lVert(w, x, y, z)\rVert = \sqrt{w^2 + x^2 + y^2 + z^2}$$.

The difference between perturbing the unit axis and rotation angle used to define the versor, and interpolating between two versors, is that the perturbation method allows you to limit the maximum angular displacement caused by the perturbation. If you only rotate the unit axis vector by angle $$\alpha$$, that is the maximum. Similarly for only modifying the rotation angle by $$\psi$$ around the axis vector. The combined limit is more complicated, but I haven't bothered to calculate it, since $$\alpha$$ and $$\psi$$ tend to be so small (less than one degree in magnitude) in practical applications that you can just assume it is less than $$\lvert\alpha\rvert+\lvert\psi\rvert$$).

In particular, to obtain a better numerical fit via trial-and-error or random walking a versor, the perturbation method should do fine. (Usually the problem is that the phase space is full of local minima: you can easily find one, but you want to find the deepest one instead. So, you need to use a "search grid" with preferably at most one local minima per grid cell, and search each grid cell separately.)

How to perturb the axis by angle $$\alpha$$, then? Pick an uniform random unit vector perpendicular to the axis, and rotate the axis vector around that by angle $$\alpha$$, of course.

• Thank you. This helps me get a clearer understanding of how to do the rejection sampling in the proper way. – Amir Nov 10 '18 at 22:04
• As Steven Stadnicki pointed out, the true uniform distribution on the set of all rotation does not yield the uniform distribution on rotation angle. Rather, you should compute the inverse function of $\beta\mapsto\beta-\sin\beta$, as explained in my answer. – Junekey Jeon Nov 11 '18 at 12:15