  This article presents a new primal-dual weak Galerkin method for second order
elliptic equations in non-divergence form. The new method is devised as a
constrained $L^p$-optimization problem with constraints that mimic the second
order elliptic equation by using the discrete weak Hessian locally on each
element. An equivalent min-max characterization is derived to show the
existence and uniqueness of the numerical solution. Optimal order error
estimates are established for the numerical solution under the discrete
$W^{2,p}$ norm, as well as the standard $W^{1,p}$ and $L^p$ norms. An
equivalent characterization of the optimization problem in term of a system of
fixed-point equations via the proximity operator is presented. An iterative
algorithm is designed based on the fixed-point equations to solve the
optimization problems. Implementation of the iterative algorithm is studied and
convergence of the iterative algorithm is established. Numerical experiments
for both smooth and non-smooth coefficients problems are presented to verify
the theoretical findings.
