Bi-factor analysis is a form of confirmatory factor analysis widely used in psychological and educational measurement. The use of a bi-factor model requires specifying an explicit bi-factor structure on the relationship between the observed variables and the group factors. In practice, the bi-factor structure is sometimes unknown, in which case, an exploratory form of bi-factor analysis is needed. Unfortunately, there are few methods for exploratory bi-factor analysis, with the exception of a rotation-based method proposed in Jennrich and Bentler ([2011, Psychometrika 76, pp. 537–549], [2012, Psychometrika 77, pp. 442–454]). However, the rotation method does not yield an exact bi-factor loading structure, even after hard thresholding. In this article, we propose a constraint-based optimization method that learns an exact bi-factor loading structure from data, overcoming the issue with the rotation-based method. The key to the proposed method is a mathematical characterization of the bi-factor loading structure as a set of equality constraints, which allows us to formulate the exploratory bi-factor analysis problem as a constrained optimization problem in a continuous domain and solve the optimization problem with an augmented Lagrangian method. The power of the proposed method is shown via simulation studies and a real data example.