We mathematically defined a higher-order GSVD for two or more large-scale matrices with different row dimensions and the same column dimension. We proved that our new HO GSVD extends to higher orders almost all of the mathematical properties of the GSVD: The eigenvalues of S are always greater than or equal to one, and an eigenvalue of one corresponds to a right basis vector of equal significance in all matrices, and to a left basis vector in each matrix CX-4945 factorization that is orthogonal to all other left basis vectors in that factorization. We therefore mathematically defined, in analogy with the GSVD, the common HO GSVD subspace of the N§2 matrices to be the subspace spanned by the right basis vectors that correspond to the eigenvalues of S that equal one. The only property that does not extend to higher orders in general is the complete column-wise orthogonality of the normalized left basis vectors in each factorization. Recent research showed that several higher-order generalizations are possible for a given matrix decomposition, each preserving some but not all of the properties of the matrix decomposition. The HO GSVD has the interesting property of preserving the exactness and diagonality of the matrix GSVD and, in special cases, also partial or even complete column-wise orthogonality. That is, all N matrix factorizations in Equation are exact, all N matrices Si are diagonal, and when one or more of the eigenvalues of S equal one, the corresponding left basis vectors in each factorization are orthogonal to all other left basis vectors in that factorization. The complete column-wise orthogonality of the matrix GSVD enables its stable computation. We showed that each of the right basis vectors that span the common HO GSVD subspace is a generalized singular vector of all pairwise GSVD factorizations of the matrices Di and Dj with equal corresponding generalized singular values for all i and j. Since the GSVD can be computed in a stable way, the common HO GSVD subspace can also be computed in a stable way by computing all pairwise GSVD factorizations of the matrices Di and Dj. That is, the common HO GSVD subspace exists also for N matrices Di that are not all of full column rank. This also means that the common HO GSVD subspace can be formulated as a solution to an optimization problem, in analogy with existing variational formulations of the GSVD. It would be ideal if our procedure reduced to the stable computation of the matrix GSVD when N~2. To achieve this ideal, we would need to find a procedure that allows a computation of the HO GSVD, not just the common HO GSVD subspace, for N matrices Di that are not all of full column rank. A formulation of the HO GSVD, not just the common HO GSVD subspace, as a solution to an optimization problem may lead to a stable numerical algorithm for computing the HO GSVD. Such a formulation may also lead to a higher-order general GaussMarkov linear statistical model. It was shown that the GSVD provides a mathematical framework for sequence-independent comparative modeling of DNA microarray data from two organisms.