You might to apply independent component analysis to your problem.  ICA has a R package call fastICA.  What happens is that you decompose X1 and X2 into two independent variables.  Then you develop your regression on the components. 

Do you mean that first PCA scores matrix is prepared and then perform regression?

Since there are M dependent variables and 2M independent variables (because of X1 and X2 matrices) then how to obtain regression equation?

No.  The ICA.   Do a google on ICA.   

You can use bayesian analysis, where the correlation between X1 and X2 is not a problem..

Thanks for your guidance.  I will look into ICA and bayesian approaches.

You can decorrelated your variables using the Cholesky decomposition of your correlation matrice Ro.

Ro=L.Lt;

X=L.U

so the variables U are incorrelted...

Perhaps You can develop a Y matrix with the sizes and ranks of these two matrices x1 and x1. 

I am thankful to the researchers for suggestions with full interest.  I feel that a complete solution is yet to arrive.  

I have no idea that how Cholesky's decomposition will work and the regression results can be interpreted specially for the situation when all the three matrices X1, X2 and Y are correlated? 

Dear Mukesh,

Since the size of three matrices X1, X2 and Y are same (NxM)  where

(N < M), and Y is related to X1 and Y is related to X2.  Suppose that

Y = A1X1  and  also  Y = A2X2  for some  A1  and A2.

Note:  A1  and  A2  must be square matrices of order NxN, by conformablity and check the orders of both sides of the above equations. 

Therefore

2Y = A1X1 + A2X2

implies

Y = (1/2)A1X1 + (1/2)A2X2

We already  develop a relationship of Y jointly with X1 and X2; i.e.

Y = A*X1 + B*X2,

where A* = (1/2)A1 and B* = (1/2)A2 .

The remainder cases, we can develop in this ways.