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.