Matrix representation for $F_4$
As I understand the question, the OP would be happy to see a description of the lowest-dimension fundamental representation of $F_4$ (and perhaps $E_6$, $E_7$, $E_8$), and is happy with the description of $G_2$ acting on the space of trace-zero octonions.
Yes - people have "bothered to write down" representations of $F_4$ and the other exceptional groups in this spirit. Chevalley, Schafer, Albert, Jacobson, Freudenthal, and Tits are the names that come to mind first.
$F_4$ acts naturally as automorphisms of the 27-dimensional exceptional Jordan algebra $J_{3,O}$ -- this is the Jordan algebra of 3 by 3 Hermitian matrices with entries in the octonions (which octonion algebra you use depends on or determines the form of $F_4$). Since $F_4$ acts as algebra automorphisms, it preserves the unit element of this algebra, and preserves the trace form as a result. It follows that $F_4$ acts on the 26-dimensional trace-zero subspace of the 27-dimensional algebra. This is quite close, in spirit, to the example of $G_2$ acting on trace-zero octonions.
Also, $F_4$ acts on this 26-dimensional space, preserving the nondegenerate symmetric trace form: $$(X,Y) \rightarrow Tr(X \cdot Y).$$
$E_6$ also acts on the 27-dimensional Jordan algebra above, but not as algebra automorphisms. Instead, $E_6$ can be viewed as the linear automorphisms of this 27-dimensional space that preserve the cubic norm form (the "determinant" of a 3 by 3 Hermitian octonionic matrix). I believe this goes back to Chevalley and Schafer about 60 years ago.
$E_7$ acts naturally on a 56-dimensional space, studied by Freudenthal. This is the space of two-by-two matrices, with diagonal entries in the base field, and off-diagonal entries in the exceptional Jordan algebra mentioned above: $2+27+27 = 56$. $E_7$ can be viewed as the group of linear automorphisms of this 56-dimensional space preserving a quartic form, I believe.
The smallest irreducible representation of $E_8$ is the adjoint representation of $E_8$ on its own Lie algebra -- so you have to construct $E_8$ to represent it, in a sense.
A nice recent survey of related topics, and a source of other references is Baez's survey on the octonions.
While there have been many people who have done this, the first person to do so was Élie Cartan, who wrote down a basis for the matrix algebra ${\frak{f}}_4\subset{\frak{so}}(26)$ in his 1894 thesis Sur la structure des groupes de transformations finis et continus. See pp 144--147 for the explicit formulae over $\mathbb{C}$ or $\mathbb{R}$ (for the split form). In a later paper, Les groupes réels simples finis et continus (1914), he explicitly exhibited (conjugate linear) involutions on $\mathbb{C}^{26}$ whose fixed subgroups would yield the three real forms of $F_4$, see pp. 343--352.
He also gives explicit formulae for $E_6$, $E_7$ and $E_8$ in their lowest dimensional representations in the above two papers.
Apropos the OP's comment My main question is really "what does the matrix look like?" What do you have to impose on the entries of the so(26) matrix to be in the F4 subalgebra?: It seems to me that you are asking how the linear transformations in the various algebras are characterized. So, for example, one might say that ${\frak{so}}(n)$ is characterized by being the skewsymmetric $n$-by-$n$ matrices or, alternatively, you might say that these are the linear transformations that, to first order, preserve the standard (positive definite) inner product on $\mathbb{R}^n$.
Along these lines, Cartan again gives the answers for $F_4$, $E_6$, and $E_7$ in his 1894 thesis: He shows that ${\frak{f}}_4\subset{\frak{so}}(26)$ is characterized as the Lie algebra of the stabilizer of a quadratic form and a certain projective cone of dimension $16$ in $26$ dimensions. (Later, he realized that it could be characterized as the stabilizer of a homogeneous cubic polynomial in $26$ variables.) He also described ${\frak{e}}_6\subset{\frak{sl}}(27)$ as the stabilizer of a homogeneous cubic polynomial in $27$ variables, and ${\frak{e}}_7\subset{\frak{sp}}(28)$ as the stabilizer of a homogeneous quartic polynomial in $56$ variables (and a symplectic form, though, for the Lie algebra, this is not necessary).
Finally, although Cartan does not seem to have noticed this, it turns out that ${\frak{e}}_8\subset{\frak{so}}(248)$ can be characterized in ${\frak{gl}}\bigl({\frak{e}}_8\bigr)$ as the stabilizer of the Cartan $3$-form of the algebra ${\frak{e}}_8$ itself.
Here I quote some GAP code for generating F4 Lie algebra as matrices 27 x 27. It is possible to obtain them in dimension 26 but then you need to use Sqrt(2) or Sqrt(3). In dimension 27 it is more nice.
Matrices L1..L7 are 8x8 matrices of left multiplication by imaginary unit octonions e1..e7. R1..R7 are right multiplications by unit octonions. If you have troubles to obtain such please let me know.
I have also obtained Lie algebra of E6 as complex matrices 27x27 and Lie algebra E7 as quaternion matrices 28x28. I tried to obtain some nice way E8 Lie algebra in shape 31*8 dim matrices but no luck. It was in 2008 when I worked on these scripts. I have used Freudenthal, Tits, Vinberg papers when doing this.
Regards, Marek
v18 := BlockMatrix([[1,1,-R1], [2,2,-L1], [3,3,L1+R1]], 4,4);
v28 := BlockMatrix([[1,1,-R2], [2,2,-L2], [3,3,L2+R2]], 4,4);
v38 := BlockMatrix([[1,1,-R3], [2,2,-L3], [3,3,L3+R3]], 4,4);
v48 := BlockMatrix([[1,1,-R4], [2,2,-L4], [3,3,L4+R4]], 4,4);
v58 := BlockMatrix([[1,1,-R5], [2,2,-L5], [3,3,L5+R5]], 4,4);
v68 := BlockMatrix([[1,1,-R6], [2,2,-L6], [3,3,L6+R6]], 4,4);
v78 := BlockMatrix([[1,1,-R7], [2,2,-L7], [3,3,L7+R7]], 4,4);
S:= DiagonalMat([1,-1,-1,-1,-1,-1,-1,-1]);
Build now elliptic version of F4. The name is from that exp(t*[[0,-1], [2,0]]) is ellipse.
Really exp(t*[[0,-1], [2,0]]/Sqrt(2)) is ellipse.
p1 := [[0, -1, 1]];
p2 := [[-1, 1, 0]];
n := NullMat(8,8);
v:= n+p1;
vt := -2*TransposedMat(v);
This a0 really corresponds to ad([[0,-1,0],[1,0,0],[0,0,0]]) derivation on h3O.
a0 := BlockMatrix([[1,2,S], [2,1,-S], [3,4, v], [4,3, vt]], 4,4);
w := n+p2;
wt := -2*TransposedMat(w);
b0 := BlockMatrix([[2,3,S], [3,2,-S], [1,4, w], [4,1, wt]], 4,4);
mats3 := [ v18, v28, v38, v48, v58,v68,v78, a0, b0];;
ms:=List(mats3, x->x{[1..27]}{[1..27]});
f4_e := LieAlgebra( Rationals, ms);
2016-05-24 I add following matrices more nice for human eye, when I learned how to entered them in MathJax. There are seven matrices generating $so_8$ (triality on Lie algebra). $$v_{k,8}=\pmatrix{ -R & & & & & \\ & -L & & & & \\ & & L+R & & & \\ & & & 0 & & \\ & & & & 0 & \\ & & & & & 0 } $$ In this notation $L$ and $R$ are matrices of left and right multiplication by imaginary base octonions. We add two more matrices to obtain $f_4$ Lie algebra. We should understand that vectors 1x8 or 8x1 have element on first position (position of 1 in octonions). Matrix S represent octonions conjugation. $$a_{0}=\pmatrix{ & S & & & & \\ -S & & & & & \\ & & 0 & & -1 & 1 \\ & & & 0 & & \\ & & 2 & & 0 & \\ & & -2 & & & 0 } $$
$$b_{0}=\pmatrix{ 0 & & & -1 & 1 & \\ & & S & & & \\ & -S & & & & \\ 2 & & & 0 & & \\ -2 & & & & 0 & \\ & & & & & 0 } $$
EDIT 2018-08-20
On request of @AccidentalFourierTransform I include multiplication table of octonions. From it one can deduce matrices $L$, $R$ of multiplication by base octonion from the left and from the right. The multiplication table which I include is obtained from quaternions by applying Cayley-Dickson construction. Exactly the same multiplication table is present on wikipedia article on octonion. I added some colors in order visualize symmetries.
This table I have created when trying to submit article to "Delta" about octonions. It is popular journal in Polish for high school pupils. I was inspired by this article by Zbigniew Marciniak. Unfortunately my article was rejected.
Regards,