Darrigol's History of Optics names Barrow (1669) and Huygens (1653) as giving semi-verbal case by case variants of the lens formula before Halley, with Molyneux being the first to publish in 1692. Halley (1693) is credited as the first who wrote the formula algebraically. Dijksterhuis in Lenses and Waves gives detailed description of Huygens's Dioptrica (which he calls Tractatus), with reproductions of his diagrams that look roughly modern, and remarks:
"In the twentieth proposition of Tractatus, Huygens extends the results for thin lenses to non-parallel rays... Huygens had to treat all cases of positive and negative lens sides separately, but the result comes down to the modern formula $\frac{1}{p'}+\frac{1}{p}=\frac{1}{f}$" [p.18].
A systematic development of modern Gaussian/paraxial optics that the formula naturally belongs to first appeared in print in Hamilton's works of 1828-33. Gauss did it earlier but his memoir was only published in 1840. Here are the relevant excerpts from Darrigol:
"Barrow, who treated optics as a mathematical playground, did not dwell on the paraxial
(Gaussian optics) that is most useful to the theory of optical instruments. He nonetheless
gave relations equivalent to the various cases of the modern formula
$$\frac{n}{p'}+\frac{1}{p}=\frac{n-1}{R}$$
relating the algebraic distances $p$ and $p'$
(from the intersection of the axis with the refracting
surface) of two conjugate points of a spherical refracting sphere with the index $n$ and the
radius $R$, as well as unnecessarily complicated relations of conjugation for thin lenses. In
contrast, Huygens’s early private treatises were aimed at improving the telescopes that he was
designing with his brother. In his Dioptrica of 1653, he derived a variant of the now familiar
$$\frac{1}{p'}+\frac{1}{p}=\frac{1}{f}\text{ with }\frac{1}{f}=(n-1)\left(\frac{1}{R}-\frac{1}{R'}\right)$$
for a thin lens of index $n$, radiuses $R$ and $R'$, and focal length $f$. This formula first appeared
in print forty years later: case by case in Molyneux’s treatise, and in algebraic form in a
memoir by Edmond Halley also including thick lenses." [pp. 75-6]
"With these powerful tools, Hamilton developed the study of caustics, foci, and various
kinds of aberrations. In a memoir read at the British Association meeting of 1833, he
briefly indicated how to use the $T$ function in the case of optical instruments with axial symmetry... The $T$ function thus leads to the basic concepts and relations of paraxial optics, although Carl Friedrich Gauss was first to define them in a general manner
in a theory conceived in his youth and published belatedly in 1840. In his demonstrations,
Gauss relied on the combination of the linear transformations connecting the parameters
(inclination and transverse coordinates in a fixed plane) of incoming and outgoing rays for
each refracting surface, as is most commonly done nowadays." [pp. 266-7]
Footnotes to the quoted text give references to the original sources.
