This paper presents a feasible method for optimization with orthogonality constraints. The method uses the Cayley transform, a Crank-Nicolson-like update scheme, to preserve constraints and develop curvilinear search algorithms with lower computational complexity compared to projection and geodesic-based methods. The algorithms are tested on various problems, including maxcut, polynomial optimization, nearest correlation matrix estimation, and quadratic assignment problems. For the maxcut problem, the algorithm exactly solves a decomposition formulation for the SDP relaxation. For polynomial optimization, nearest correlation matrix estimation, and extreme eigenvalue problems, the algorithms run efficiently and return solutions no worse than state-of-the-art methods. For the quadratic assignment problem, the algorithm achieves a gap of 0.842% to the best known solution on the largest problem "tai256c" in QAPLIB in 5 minutes on a typical laptop. The paper considers optimization problems with orthogonality constraints, such as minimizing a function subject to X^T X = I, and generalized orthogonality constraints X*MX = K. These constraints are challenging due to their non-convexity and computational expense. Existing constraint-preserving algorithms often use matrix re-orthogonalization or geodesic paths, which require matrix factorizations or solving PDEs. The proposed method uses the Cayley transform to preserve constraints and is efficient for both real and complex optimization problems. The results apply to more general problems after trivial modifications. The paper highlights the efficiency and effectiveness of the proposed algorithms in solving optimization problems with orthogonality constraints.This paper presents a feasible method for optimization with orthogonality constraints. The method uses the Cayley transform, a Crank-Nicolson-like update scheme, to preserve constraints and develop curvilinear search algorithms with lower computational complexity compared to projection and geodesic-based methods. The algorithms are tested on various problems, including maxcut, polynomial optimization, nearest correlation matrix estimation, and quadratic assignment problems. For the maxcut problem, the algorithm exactly solves a decomposition formulation for the SDP relaxation. For polynomial optimization, nearest correlation matrix estimation, and extreme eigenvalue problems, the algorithms run efficiently and return solutions no worse than state-of-the-art methods. For the quadratic assignment problem, the algorithm achieves a gap of 0.842% to the best known solution on the largest problem "tai256c" in QAPLIB in 5 minutes on a typical laptop. The paper considers optimization problems with orthogonality constraints, such as minimizing a function subject to X^T X = I, and generalized orthogonality constraints X*MX = K. These constraints are challenging due to their non-convexity and computational expense. Existing constraint-preserving algorithms often use matrix re-orthogonalization or geodesic paths, which require matrix factorizations or solving PDEs. The proposed method uses the Cayley transform to preserve constraints and is efficient for both real and complex optimization problems. The results apply to more general problems after trivial modifications. The paper highlights the efficiency and effectiveness of the proposed algorithms in solving optimization problems with orthogonality constraints.