Getting the arbitrary boost matrix from a similarity transformation Note: For the following question I'm using the non-standard $(x,y,z,ct)$ notation. 
I'm wanting to represent an arbitrary boost in the $\hat{\beta}$ direction by doing a similarity transformation on the frame boosted in the $z$ direction to a frame with axes $\{{r_1,r_2,\hat{\beta}}\}$, allowing $r_1$ and $r_2$ to be any two mutually orthonormal vectors, as well as orthogonal to $\hat{\beta}$. So the boost matrix and transformation matrix are, respectively, 
$$L=\pmatrix{1&0&0&0\\0&1&0&0\\0&0&\gamma&-\gamma\beta\\0&0&-\gamma\beta&\gamma}$$
$$R = \pmatrix{x\cdot r_1&y\cdot r_1&z\cdot r_1&0\\x\cdot r_2&y\cdot r_2&z\cdot r_2&0\\x\cdot\hat{\beta}&y\cdot\hat{\beta}&z\cdot\hat{\beta}&0\\0&0&0&1} = \pmatrix{r_{1x}&r_{1y}&r_{1z}&0\\r_{2x}&r_{2y}&r_{2z}&0\\\beta_{x}/\beta&\beta_{y}/\beta&\beta_{z}/\beta&0\\0&0&0&1}$$
Carrying out the similarity transformation,
$$\tilde{L}=RLR^T=\pmatrix{r_{1x}^2+r_{1y}^2+\gamma r_{1z}^2&r_{1x}r_{2x}+r_{1y}r_{2y}+\gamma r_{1z}r_{2z}&(r_{1x}\beta_{x}+r_{1y}\beta_{y}+\gamma r_{1z}\beta_{z})/\beta&-\gamma r_{1z}\beta\\r_{1x}r_{2x}+r_{1y}r_{2y}+\gamma r_{1z}r_{2z}&r_{2x}^2+r_{2y}^2+\gamma r_{2z}^2&(r_{2x}\beta_{x}+r_{2y}\beta_{y}+\gamma r_{2z}\beta_{z})/\beta&-\gamma r_{2z}\beta\\(r_{1x}\beta_{x}+r_{1y}\beta_{y}+\gamma r_{1z}\beta_{z})/\beta&(r_{2x}\beta_{x}+r_{2y}\beta_{y}+\gamma r_{2z}\beta_{z})/\beta&(\beta_x^2+\beta_y^2+\gamma\beta_z^2)/\beta^2&-\gamma\beta_z\\-\gamma r_{1z}\beta&-\gamma r_{2z}\beta&-\gamma\beta_z&\gamma}$$
Utilizing mutual orthogonality to eliminate all components of the $r$'s except $r_{1z}$ and $r_{2z}$ results in 
$$\tilde{L}=\pmatrix{1+r_{1z}^2(\gamma-1)&r_{1z}r_{2z}(\gamma-1)&r_{1z}\beta_z(\gamma-1)/\beta&-\gamma\beta r_{1z}\\r_{1z}r_{2z}(\gamma-1)&1+r_{2z}^2(\gamma-1)&r_{2z}\beta_z(\gamma-1)/\beta&-\gamma\beta r_{2z}\\r_{1z}\beta_z(\gamma-1)/\beta&r_{2z}\beta_z(\gamma-1)/\beta&1+\beta_z^2(\gamma-1)/\beta^2&-\gamma\beta_z\\-\gamma\beta r_{1z}&-\gamma\beta r_{2z}&-\gamma\beta_z&\gamma}$$
And the end desired result from the exercise is
$$\tilde{L} = RLR^T = \pmatrix{1+\frac{\beta_{x}^2(\gamma-1)}{\beta^2}&\frac{\beta_{x}\beta_{y}(\gamma-1)}{\beta^2}&\frac{\beta_{x}\beta_{z}(\gamma-1)}{\beta^2}&-\beta_{x}\gamma\\\frac{\beta_{x}\beta_{y}(\gamma-1)}{\beta^2}&1+\frac{\beta_{y}^2(\gamma-1)}{\beta^2}&\frac{\beta_{y}\beta_{z}(\gamma-1)}{\beta^2}&-\beta_y\gamma\\\frac{\beta_{x}\beta_{z}(\gamma-1)}{\beta^2}&\frac{\beta_{y}\beta_{z}(\gamma-1)}{\beta^2}&1+\frac{\beta_{z}^2(\gamma-1)}{\beta^2}&-\beta_z\gamma\\-\beta_x\gamma&-\beta_y\gamma&-\beta_z\gamma&\gamma}$$
My confusion: the exercise states $r_1$ and $r_2$ to be any two orthonormal vectors, but I specifically need them to have their $z$-components be  $r_{1z}=\beta_x/\beta$ and $r_{2z}=\beta_y/\beta$ to get agreement with the matrix above. How am I going about this in the wrong way so that the result is dependent on these $z$-components of $r_1$ and $r_2$? Perhaps I'm supposed to utilize invariance of the determinant under similarity transformations, $|\tilde{L}|=|L|$, to eliminate one more of those $z$-components? That just seems like an algebraic nightmare, and with my matrix being so close in form I don't think it's the correct approach. 
 A: The general Lorentz transformation matrix is:
$$L=\begin{bmatrix}
  \gamma & -\gamma\,\vec{\beta}^T \\
  -\gamma\,\vec{\beta} & I_3+\frac{\gamma-1}{\vec{v}\,\cdot \vec{v}}\,\vec{\beta}\,\vec{\beta}^T \\
\end{bmatrix}$$
if you rotate the Lorentz matrix then:
$\vec{\beta}\mapsto R\,\vec{\beta}$ 
where R is $3\times 3$ orthogonal  Rotation Matrix
$R^T=\left[\vec{r_1},\vec{r_2},\vec{r_3}\right]$
and $R^\,R^T=I_3$
because $R\,\vec{\beta}$ in your case is :
$$\beta_z\,\begin{bmatrix}
   r_{1z} \\
   r_{2z} \\
   r_{3z} \\
 \end{bmatrix}$$
only the z-components  of the vectors that create the rotation matrix R are involved    
A: I confess I am not completely comfortable with this problem. (Witness my first attempt at an answer.)  
Pressing ahead regardless, I think that $\vec{r_1}$ and $\vec{r_2}$ are constrained by the requirement that the Lorentz transformation $\tilde{L}$ be a pure boost; that is, that spatial vectors perpendicular to the boost direction $\vec{\beta}/\beta$  are left unchanged by $\tilde{L}$.  (These vectors form a 2-dimensional eigenspace with eigenvalue 1.)
Specifically, solving $\tilde{L} \vec{r_1} = \vec{r_1}$ and $\tilde{L} \vec{r_2} = \vec{r_2}$ (or components thereof) yields 2 conditions which are satisfied by choosing $r_{1z}=\beta_x/\beta$ and $r_{2z}=\beta_y/\beta$.  
For example, the $x$-component of $\tilde{L} \vec{r_1} = \vec{r_1}$ simplifies to:
$$ r_{1x} r_{1z} + r_{1y} r_{2z} + r_{1z} \beta_z/\beta = 0 $$
With the substitutions, this equation becomes an inner product of $r_1$ and $\beta$, which vanishes by the orthogonality of $r_1$ and $\beta$.  
The other components yield the same or a similar equation which is satisfied by the same substitutions.
UPDATE:
Things are clearer now.  My answer above was problematic in that it put some additional constraints on $\hat{r_1}$ and $\hat{r_2}$ beyond their being orthonormal to each other and $\hat{\beta}$, while the problem stated that any orthonormal pair should work.  In fact, any such pair $\hat{r_1}$ and $\hat{r_2}$ yields a symmetric $\tilde{L}$, and hence a pure boost, so why should additional constraints be required?
The answer is that additional constraints are not required, if one chooses the matrix $R$ to rotate $\hat{z}$ into $\hat{\beta}$ ($R \hat{z}=\hat{\beta}$). Specifically, $R$ should have been the transpose of the one you used.  With this change, the calculation you performed gives the expected result, with no additional constraints on $\hat{r_1}$ and $\hat{r_2}$ required.  
(Conversely, with the version of $R$ you employed, it was necessary that $r_{1z}=\beta_x/\beta$ and $r_{2z}=\beta_y/\beta$ so that $\hat{z}$ was rotated into $\hat{\beta}$.)
