The Pfaffian is an invariant of matrices, but not of underlying linear operators. Two skew-symmetric matrices that are similar can have different Pfaffians: for example, $$\mathrm{Pf}\begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} = 1, \; \; \mathrm{Pf} \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} = -1.$$ In other words, the Pfaffian depends not only on the linear operator, but also on a choice of basis - so there can be no basis-free definition.


The Pfaffian only makes sense for skew-symmetric matrices in even dimensions, so there's no hope of generalizing it as far as an arbitrary linear operator.

One coordinate-free description of the input to the Pfaffian is that it is an element of the orthogonal Lie algebra $\mathfrak{so}(2n)$; in this incarnation the Pfaffian appears as an invariant polynomial on $\mathfrak{so}(2n)$ (so it is invariant under orthogonal, rather than arbitrary, changes of coordinates). It is in fact the polynomial which corresponds via Chern-Weil theory to the Euler class; this is more or less the content of the Chern-Gauss-Bonnet theorem. The Pfaffian squaring to the determinant corresponds via Chern-Weil theory to the Euler class squaring to the top Pontryagin class.

It should be possible to give a conceptual proof that the Pfaffian squares to the determinant using the wedge product definition. As in user357105's comment to his answer, the setting to work in is a finite-dimensional real (for simplicity) inner product space $V$; here you can identify skew-symmetric linear operators $T : V \to V$ with elements of $\wedge^2 V$ (both of which can in turn be identified with the Lie algebra $\mathfrak{so}(V)$), and you use this identification to relate the determinant and the Pfaffian. I haven't worked through the details, though.