I'm assuming Strang is allowing complex eigenvectors. The statement "$A$ has orthogonal eigenvectors" is a bit imprecise; a better statement would be "You can choose a basis of pairwise orthogonal eigenvectors". If you take his statement literally then it is correct but trivial, since a large enough normal matrix ($AA^*=A^* A$) will have some pair of orthogonal eigenvectors.
The eigenspaces for distinct eigenvalues are orthogonal, but that's different from saying there exists an orthonormal basis consisting of eigenvectors. The eigenspaces might not sum to the whole space.
Consider the usual dot product
$$
(x_{1}, x_{2}, x_{3}, x_{4})
\cdot
(y_{1}, y_{2}, y_{3}, y_{4})
=
x_{1} y_{1} + x_{2} y_{2} + x_{3} y_{3} + x_{4} y_{4}.
$$
Then it is easy to see that $A$ is symmetric iff for all $x, y$
$$
x A \cdot y = x \cdot y A.
$$
(just take the $x = e_{i}, y = e_{j}$, where $e_{i}$ is the vector which is all zero except for a $1$ in the $i$-th position).
Now prove that if $x$ is an eigenvector with respect to $\lambda$, $y$ is an eigenvector with respect to $\mu$, and $\lambda \ne \mu$, then $x, y$ are orthogonal:
$$
\lambda (x \cdot y)
=
(\lambda x) \cdot y
=
x A \cdot y
=
x \cdot y A
=
x \cdot (\mu y)
=
\mu (x \cdot y)
$$
and $\lambda \ne \mu$ implies $x \cdot y = 0$.
Therefore in your case the missing eigenvector $v$ relative to the eigenvalue $2$ must be orthogonal to both $(1, 0, 0, −1)$ and $(0, 1, 1, 0)$. It is easy to check that this means
$$
v = (a, b, -b, a)
=
a (1, 0, 0, 1) + b (0, 1, -1, 0),
$$
for some $a, b$.
Best Answer
As I think you are asking for intuition regarding "The second derivative in a specific direction represented by a unit vector d is given by $d’Hd$”, let me correlate it in two ways with the normal way we think about derivatives. I’ll use two dimensions to illustrate in both cases. Let the unit vector $\bar{d}$ be $(n_1,n_2)$ in the standard basis and let $\bar{x}$ represent the point (x,y).
For the shorter explanation, consider the function value $f(\bar{x}+ds \bar{d})$ at a small distance $ds$ from $\bar{x}$ along $\bar{d}$ as a Taylor expansion. Let $h=n_1ds$ and $k=n_2ds$ denote the corresponding increments along the x and y directions.
$$f(\bar{x}+ds \bar{d})=f(x,y) + hf_x+kf_y + \frac{1}{2}(h^2f_{xx}+ 2hkf_{xy}+ k^2f_{yy}) + \mbox{h.o.t.}$$
$$=f(x,y) + ds(n_1f_x+n_2f_y) + \frac{1}{2}ds^2(n_1^2f_{xx}+ 2n_1n_2f_{xy}+ n_2^2f_{yy}) + \mbox{h.o.t.}$$
$$=f(x,y) + ds (\nabla f \cdot \bar{d} )+ \frac{1}{2}ds^2 (\bar{d}’H \bar{d} )+ \mbox{h.o.t.}$$
That is, $\nabla f \cdot \bar{d}$ plays the role of the first derivative and $\bar{d}’H \bar{d}$ plays the role of the second derivative along the direction $\bar{d}$.
The second explanation is using the same idea but depending on your bent of mind, might be more intuitive. Proceeding as in finite differences, where $f_x$ is approximated by $\frac{f(x+\Delta x)-f(x)}{\Delta x}$ with the approximation becoming exact as $\Delta x \rightarrow 0$. Then the second derivative $f_{xx}$ is likewise approximated by $$\frac{ f_x(x+\frac{\Delta x}{2}) - f_x(x -\frac{\Delta x}{2}) }{\Delta x}$$
$$~ \frac{ f( x + \Delta x) -2f(x) + f( x - \Delta x) }{\Delta x^2}$$ Now, apply that one dimensional second derivative idea along the direction $\bar{d}$ to see that, ignoring higher order terms for now, the second derivative is
$$ \frac{ f( x + h, y+ k) -2f(x) + f( x - h, y-k }{ h^2 + k^2}$$
Using 2 dimensional Taylor expansions for $f( x + h, y+ k)$ and $ f( x - h, y-k )$ (write it out)
and using $h=n_1ds$ and $k=n_2ds$, we see that the second derivative approximation is given by
$$ds^2 \frac{ n_1^2f_{xx}+ 2n_1n_2f_{xy}+ n_2^2f_{yy} }{ ds^2} = ds^2 \frac{ \bar{d}’H \bar{d} }{ ds^2} = \bar{d}’H \bar{d} $$ with the higher order terms vanishing as you take $ds$ to zero.
I would have liked to expand some of the steps more, but MathJax on a phone is rather painful. I hope one of these explanations felt intuitive to you. Please leave a comment if more clarification is needed.