Abstract <p>The Moore–Penrose inverse is the most popular type of matrix generalized inverses that has many applications both in matrix theory and numerical linear algebra. The Moore–Penrose inverse of a matrix can be found via singular value decomposition (abbreviated SVD) of this matrix. In this regard, there exist the most effective algorithm which consists of two phases. At the first phase an initial matrix is reduced to upper bidiagonal form (the Golub–Kahan bidiagonalization algorithm). The second phase is known as the Golub–Reinsch algorithm. This is an iterative procedure that generates a sequence of bidiagonal matrices converging to a diagonal form. Acting in this way, we obtain an iterative approximation to the SVD of bidiagonal matrix. In this paper, we develop a method that can be considered as an alternative to the Golub–Reinsch iterative procedure. By implementing an approach proposed in the work, the following two main results have been achieved. First, we obtain explicit expressions for the entries of the Moore–Penrose inverse of bidigonal matrices. Secondly, based on the closed form formulas, we get a finite recursive numerical algorithm of optimal computational complexity. Thus, in certain cases, we can compute the Moore–Penrose inverse of bidiagonal matrices without using the SVD.</p>