Although researchers desire to evaluate system reliability accurately and efficiently over the years, little progress has been made on system reliability analysis. Up to now, bound methods for system reliability prediction have been dominant. However, two primary challenges are as follows: (1) Most numerical methods cannot effectively evaluate the probabilities of the second (or higher)–order joint failure events with high efficiency and accuracy, which are needed for system reliability evaluation and (2) there is no unique system reliability approximation formula, which can be evaluated efficiently with commonly used reliability methods. Thus, this paper proposes the complementary intersection (CI) event, which enables us to develop the complementary intersection method (CIM) for system reliability analysis. The CIM expresses the system reliability in terms of the probabilities of the CI events and allows the use of commonly used reliability methods for evaluating the probabilities of the second–order (or higher) joint failure events efficiently. To facilitate system reliability analysis for large-scale systems, the CI-matrix can be built to store the probabilities of the first- and second-order CI events. In this paper, three different numerical solvers for reliability analysis will be used to construct the CI-matrix numerically: first-order reliability method, second-order reliability method, and eigenvector dimension reduction (EDR) method. Three examples will be employed to demonstrate that the CIM with the EDR method outperforms other methods for system reliability analysis in terms of efficiency and accuracy.