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 the specification of 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 to find the bi-factor structure. Unfortunately, there are few methods for exploratory bi-factor analysis, with the exception of a rotation-based method proposed in Jennrich and Bentler (2011, 2012). However, this method only finds approximate bi-factor structures, as it does not yield an exact bi-factor loading structure, even after applying hard thresholding. In this paper, we propose a constraint-based optimisation 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 characterisation 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 optimisation problem in a continuous domain and solve the optimisation problem with an augmented Lagrangian method. The power of the proposed method is shown via simulation studies and a real data example. Extending the proposed method to exploratory hierarchical factor analysis is also discussed. The codes are available on ``https://anonymous.4open.science/r/Bifactor-ALM-method-757D".