  This article proposes a new numerical algorithm for second order elliptic
equations in non-divergence form. The new method is based on a discrete weak
Hessian operator locally constructed by following the weak Galerkin strategy.
The numerical solution is characterized as a minimization of a non-negative
quadratic functional with constraints that mimic the second order elliptic
equation by using the discrete weak Hessian. The resulting Euler-Lagrange
equation offers a symmetric finite element scheme involving both the primal and
a dual variable known as the Lagrange multiplier, and thus the name of
primal-dual weak Galerkin finite element method. Optimal order error estimates
are derived for the finite element approximations in a discrete $H^2$-norm, as
well as the usual $H^1$- and $L^2$-norms. Some numerical results are presented
for smooth and non-smooth coefficients on convex and non-convex domains.
