Partial　least　squares　regression　(PLSR)　has　been　a　popular　technique　to　explore　the　linear　relationship　between　two　data　sets.　However,　all　existing　approaches　often　optimize　a　PLSR　model　in　Euclidean　space　and　take　a　successive　strategy　to　calculate　all　the　factors　one　by　one　for　keeping　the　mutually　orthogonal　PLSR　factors.　Thus,　a　suboptimal　solution　is　often　generated.　To　overcome　the　shortcoming,　this　paper　takes　statistically　inspired　modification　of　PLSR　(SIMPLSR)　as　a　representative　of　PLSR,　proposes　a　novel　approach　to　transform　SIMPLSR　into　optimization　problems　on　Riemannian　manifolds,　and　develops　corresponding　optimization　algorithms.　These　algorithms　can　calculate　all　the　PLSR　factors　simultaneously　to　avoid　any　suboptimal　solutions.　Moreover,　we　propose　sparse　SIMPLSR　on　Riemannian　manifolds,　which　is　simple　and　intuitive.　A　number　of　experiments　on　classification　problems　have　demonstrated　that　the　proposed　models　and　algorithms　can　get　lower　classification　error　rates　compared　with　other　linear　regression　methods　in　Euclidean　space.　We　have　made　the　experimental　code　public　at　https://github.com/Haoran2014.