Motivated by advances in PDE-constrained optimization, a fundamentally new auto- nomic closure for large eddy simulations (LES) is presented that implements an optimiza- tion formulation for the subgrid-scale stresses instead of using a predefined turbulence model. The autonomic closure approach is based on the most general dimensionally- consistent polynomial expansion of the local subgrid-scale stress tensor in terms of the resolved scale primitive variables and their products at all spatial locations and times. In so doing, the closure approach inherently addresses nonlinear, nonlocal, and nonequilibrium turbulence effects without introducing any tuning parameters. The expansion introduces a large set of coefficients that can be determined by solving an inverse problem that mini- mizes error relative to known subgrid stresses at a test filter scale. The resulting optimized coefficients are then projected to the LES scale by invoking scale similarity in the inertial range and applying appropriate renormalizations. This new closure approach avoids the need to specify a subgrid-scale model, and instead allows the optimization procedure to de- termine the best local relation between subgrid stresses and resolved-scale variables. Here we present the most general formulation of this new autonomic approach, and also present an inverse approach for determining the optimal coefficients. We then explore truncation, regularization, and sampling within the inverse formulation. Finally, we present results from a priori tests of the autonomic closure approach using data from direct numerical simulations of homogeneous isotropic and sheared turbulence. Even for the simplest 2nd order truncation of the fundamental polynomial expansion, substantial improvements over the Dynamic Smagorinsky model are found from this new autonomic closure approach.