Simultaneous Diagonalization of two bilinear forms

As I commented, a more difficult version is at Congruence and diagonalizations

This one is easier, the matrix $C$ has three distinct eigenvalues that are integers. $$ A = \left( \begin{array}{rrr} 1 & \frac{1}{2} & 0 \\ \frac{1}{2} & 1 & -\frac{1}{2} \\ 0 & -\frac{1}{2} & 1 \end{array} \right) $$

$$ B = \left( \begin{array}{rrr} 0 & -2 & 4 \\ -2 & 1 & 2 \\ 4 & 2 & 0 \end{array} \right) $$

$$ A^{-1} = \left( \begin{array}{rrr} \frac{3}{2} & -1 & -\frac{1}{2} \\ -1 & 2 & 1 \\ -\frac{1}{2} & 1 & \frac{3}{2} \end{array} \right) $$

$$ C = A^{-1}B = \left( \begin{array}{rrr} 0 & -5 & 4 \\ 0 & 6 & 0 \\ 4 & 5 & 0 \end{array} \right) $$

The theorem from Horn and Johnson (first edition hardcover was 1985, paperback 1990) is that we can continue if and only if $C$ is diagonalizable in that $R^{-1} C R =D$ is diagonal.

They were extremely careful: the eigenvalues of $C$ are $6,4,-4$ and we can make the matrix $R$ with columns as eigenvectors with $$ R = \left( \begin{array}{rrr} 1 & 1 & -1 \\ -2 & 0 & 0 \\ -1 & 1 & 1 \end{array} \right). $$ Confirm $$ CR = \left( \begin{array}{rrr} 6 & 4 & 4 \\ -12 & 0 & 0 \\ -6 & 4 & -4 \end{array} \right). $$ It follows that $R^{-1}CR$ is the diagonal matrix with entries $6,4,-4.$ $$ D = \left( \begin{array}{rrr} 6 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & -4 \end{array} \right). $$

Finally, $$ R^TBR = \left( \begin{array}{rrr} 12 & 0 & 0 \\ 0 & 8 & 0 \\ 0 & 0 & -8 \end{array} \right), $$ and $$ R^TAR = \left( \begin{array}{rrr} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{array} \right). $$ The reason this works is that we have arranged $R^T AR D = R^T BR.$

Let's see, this is in the first edition of Horn and Johnson, table 4.5.15T on page 229, then detail for case II on pages 231-232. The technique gives the full problem when $C$ has all eigenvalues distinct, as in this problem. With repeat eigenvalues, one needs to continue working, that is half of page 232. Indeed, when this happens, we are guaranteed to have square diagonal blocks in both the revised $A$ and $B$ matrices, such that each $B$ block is just a scalar multiple of each corresponding $A$ block. The extra work is then to diagonalize the $A$ block, the same thing will work on the $B$ block. That is exactly what is shown at Congruence and diagonalizations

Oh, you wanted the identity matrix. Do another step with $R_2 = I/ \sqrt 2$ on right and left, the result is just to halve each diagonal matrix. The final overall matrix is my $R/\sqrt 2.$

Tuesday: to continue, take $Q = R/\sqrt 2.$ Then $Q^{-1} = R^{-1} \sqrt 2$ $$ Q^{-1} = R^{-1} \sqrt 2 = \; \; \left( \begin{array}{rrr} 0 & \frac{-1}{\sqrt 2} & 0 \\ \frac{1}{\sqrt 2} & 0 & \frac{1}{\sqrt 2} \\ \frac{-1}{\sqrt 2} & \frac{-1}{\sqrt 2} & \frac{1}{\sqrt 2} \end{array} \right). $$

The rows of $Q^{-1}$ give us $$ \left(\frac{-1}{\sqrt 2}y \right)^2 + \left(\frac{1}{\sqrt 2} x + \frac{1}{\sqrt 2} z \right)^2 + \left(\frac{-1}{\sqrt 2} x +\frac{-1}{\sqrt 2}y + \frac{1}{\sqrt 2} z \right)^2 = x^2 + y^2 + z^2 -yz +xy $$

$$ 6\left(\frac{-1}{\sqrt 2}y \right)^2 + 4 \left(\frac{1}{\sqrt 2} x + \frac{1}{\sqrt 2} z \right)^2 -4 \left(\frac{-1}{\sqrt 2} x +\frac{-1}{\sqrt 2}y + \frac{1}{\sqrt 2} z \right)^2 = y^2 + 4yz + 8zx - 4 xy $$

(This is too long for a comment.)

The OP says in a comment to another answer here that they can only use "Lagrange diagonalization". I have never heard of this term before. From the very limited sources I could find on this site and on the internet, this seems to be just a recursive method that performs repeatedly rank-1 updates to the quadratic form, so that the rank of the updated form is strictly decreasing, until nothing remains.

E.g. the matrix representation of your $f$ is given by $$ A = \pmatrix{1&\frac12&0\\ \frac12&1&-\frac12\\ 0&-\frac12&1}. $$ The first column of $A$ is $u=(1,\frac12,0)^T$, whose first entry is $u_1=1$. Hence we split $A$ into the sum of $A$ and $\frac{1}{u_1}uu^T$ and $A-\frac{1}{u_1}uu^T$, and we continue in this manner: \begin{align} A&=1\pmatrix{1\\ \frac12\\ 0}\pmatrix{1&\frac12&0} - \pmatrix{0&0&0\\ 0&1&2\\ 0&2&4}\\ &=1\pmatrix{1\\ \frac12\\ 0}\pmatrix{1&\frac12&0} -1\pmatrix{0\\ 1\\ 2}\pmatrix{0&1&2} + 0 \end{align} and therefore $$ A=\underbrace{\pmatrix{1&0&\ast\\ \frac12&1&\ast\\ 0&2&\ast}}_{P^T} \pmatrix{1\\ &1\\ &&0} \underbrace{\pmatrix{1&\frac12&0\\ 0&1&2\\ \ast&\ast&\ast}}_P, $$ where the unspecified entries (marked by the asterisks) are arbitrary as long as $P$ is invertible.

As each iteration of this procedure depends entirely on the first column of the remaining term, apparently it is not possible to use this procedure to simultaneously diagonalise two symmetric matrices if their first columns are linearly independent, and this is your case because the matrix for $g$ is given by $$ B=\pmatrix{0&-2&4\\ -2&1&2\\ 4&2&0}. $$ But that doesn't mean simultaneous diagonalisation is impossible. It's just that Lagrange diagonalisation doesn't work. See the other answer here.

I wonder if the method the OP was attempting is simply Lagrange's method from multivariable calculus. Because a quadratic form takes its extreme values on the unit circle (or any circle) at eigenvectors, and because we can diagonalize a symmetric matrix using a basis of orthonormal eigenvectors, we look for the stationary points of $g$ with the constraint $f(x,y,z)=1$ (here $f$ will be our scalar product since, as the OP mentioned, it is positive definite and this is enough to guarantee a solution). I'll represent $f$ and $g$ by matrices $A$ and $B$ respectively to match the other answers, so we have $$A=\begin{bmatrix}1&\frac{1}{2}&0\\\frac{1}{2}&1&-\frac{1}{2}\\0&-\frac{1}{2}&1\\\end{bmatrix} \qquad B=\begin{bmatrix}0&-2&4\\-2&1&2\\4&2&0\\\end{bmatrix}$$ Lagrange's method results in the system of linear equations $$\begin{bmatrix}-\lambda&-2-\lambda/2&4\\-2-\lambda/2&1-\lambda&2+\lambda/2\\4&2+\lambda/2&-\lambda\\\end{bmatrix}\begin{bmatrix}x\\y\\z\\\end{bmatrix}=\begin{bmatrix}0\\0\\0\\\end{bmatrix}$$ where the $(i,j)$-th entry in the leftmost matrix is $b_{ij}-\lambda a_{ij}$. The system has a nontrivial solution when that matrix's determinant vanishes, which is when $\lambda=6$ or $\lambda=4$ or $\lambda=-4$. Now, we can plug in each of these values in turn to the system to get solutions, which will be eigenvectors forming a basis $\{b_1,b_2,b_3\}$ that makes both quadratic forms diagonal. I get $$b_1=\begin{bmatrix}-1\\2\\1\\\end{bmatrix},\quad b_2=\begin{bmatrix}1\\0\\1\\\end{bmatrix},\quad b_3=\begin{bmatrix}-1\\0\\1\\\end{bmatrix}$$ The final step is to normalize these vectors so that $A$ is not only diagonal, but the identity matrix. We need to do this because the information that $f(x,y,z)=1$ was lost when taking partial derivatives; we could have started out with $f(x,y,z)=\pi^3$ and got the same system of equations, and furthermore I arbitrarily chose $b_1,b_2,b_3$ to have third coordinate $1$ so there's no reason to think they'd be normalized. Just keep in mind that our scalar product is $f$ and not the typical dot product, so in fact every vector in the above basis has magnitude $\sqrt{2}$. If we put all this together, our change of basis matrix is $$R=\dfrac{1}{\sqrt{2}}\begin{bmatrix}-1&1&-1\\2&0&0\\1&1&1\\\end{bmatrix}$$ and the matrices $\widetilde{A}$, $\widetilde{B}$ of the quadratic forms in the new basis are $$\widetilde{A}=R^TAR=I \qquad \widetilde{B}=R^TBR=\begin{bmatrix}6&0&0\\0&4&0\\0&0&-4\\\end{bmatrix}$$ Two remarks: This method is described in Shilov's Linear Algebra, section 10.32. Also, I've been referring to column matrices as vectors, which if we're being pedantic is tacitly assuming that we're starting in the standard basis. However, the method works no matter what basis we start with. If it's not the standard basis, then the columns of the change of basis matrix simply represent the coordinates of the new basis vectors in terms of the starting basis, instead of the vectors themselves.