回転の表現を「標準化」する

多次元の回転について、前回予想していたことが確認できてしまいました。数学的な証明ではありませんが、実際、計算できています。

乱数によって n 次元空間に n 個のベクトルを置きます。それを一個ずつ、それ以前のベクトルのどれとも直交させ（それは内積部分を差引くことでできますが）長さを 1 に合わせることによって、直交基底 O （オー）を作ります。このものを、O = L ・ M ・ R という形に分解します。L、M、R は何れも直交行列であり、L と R は互いに転置行列（この場合逆行列）です。過去にご説明した「ODO 分解」に似ており、あの時の真ん中は対角行列でしたが、今度の M は規格化された直交行列になります。また、あの時の左右は必ずしも逆行列（転置行列）ではなかったですが、今度のは互いに転置行列です。ランダムに作った 2 から 9 次元の行列の結果として、M は次のようになります。

狙ったわけではないのですが、このように対角成分が大きいものから順に並ぶのは不思議です。

プログラムです。L も、結果は示しませんでしたが計算できています。 R は出していませんが、L を転置させるだけで良いので。

using System;
	using System.Drawing;
	using System.Windows.Forms;
	 
	public class MyForm : Form { 
		[STAThread]
		static void Main() {  
			Application.Run(new MyForm());
		} // ================================================================= 
		double[][] createOrthogonal(int n) { // === create a random orthogonal matrix !!! 
			double[][] a=new double[n][];
			Ran0 ran=new Ran0(345); 
				// === 0 to 1; random number generator; later in this file !!! 
			for (int i=0; i<n; i++) {
				a[i] = new double[n]; a[i][i] = 1; 
				for (int j=0; j<n; j++) {
					double d = ran.step();
					a[i][j] += (d - 0.5) / 1.5; } } 
			for (int i=0; i<n; i++) {
				double sum1 = 0; 
				for (int j=0; j<n; j++) {
					sum1 += a[i][j] * a[i][j]; }
				double fct = 1.0 / Math.Sqrt(sum1); 
				for (int j=0; j<n; j++) {
					a[i][j] *= fct; } 
				for (int j=i+1; j<n; j++) {
					double sum2 = 0; 
					for (int k=0; k<n; k++) {
						sum2 += a[i][k] * a[j][k]; } 
					for (int k=0; k<n; k++) {
						a[j][k] -= sum2*a[i][k]; } } } 
			return a; 
		} // ===================================================================== 
		void zeroaverageij(double[][] L, double[][] mat, int i, int j) { 
				// === average of mat[i][j] and mat[j][i] will be closer to zero !!! 
			double tt=Math.Atan2(mat[i][j]+mat[j][i],mat[i][i] - mat[j][j])/2, 
					C=Math.Cos(tt), S=Math.Sin(tt);
			for (int k=i; k<mat.Length; k++) {
				double a = mat[i][k], b = mat[j][k]; 
				double an = a * C + b * S, bn = b * C - a * S;
				mat[i][k] = an;  mat[j][k]=bn; } 
			for (int k=i; k<mat.Length; k++) {
				double a = mat[k][i], b = mat[k][j]; 
				double an = a * C + b * S, bn = b * C - a * S;
				mat[k][i] = an;  mat[k][j]=bn; } 
			for (int k=0; k<mat.Length; k++) {
				double a = L[k][i], b = L[k][j]; 
				double an = a * C + b * S, bn = b * C - a * S;
				L[k][i] = an;  L[k][j]=bn; } 
		} // =====================================================================  
		void equalizeij(double[][] L, double[][] mat, int i, int j) {
				// === distinction between mat[i][j] and mat[j][i] will be zero !!! 
			int i1 = i + 1;  
			double tt=Math.Atan2(mat[i][j] - mat[j][i],mat[i][i1] - mat[i1][i]), 
					C=Math.Cos(tt), S=Math.Sin(tt);
			for (int k=i; k<mat.Length; k++) {
				double c = mat[i1][k], d = mat[j][k];
				double cn = c * C + d * S, dn = d * C - c * S;  
				mat[i1][k] = cn; mat[j][k] = dn; } 
			for (int k=i; k<mat.Length; k++) {
				double a = mat[k][i1], b = mat[k][j];
				double an = a * C + b * S, bn = b * C - a * S;
				mat[k][i1] = an;  mat[k][j] = bn; } 
			for (int k=0; k<mat.Length; k++) { 
				double c = L[k][i1], d = L[k][j]; 
				double cn = c * C + d * S, dn = d * C - c * S; 
				L[k][i1] = cn; L[k][j] = dn; } 
		} // ===================================================================== 
		void getregular(double[][] L, double[][] mat) { 
			int i=0; 
			if (mat.Length%2 == 1) { 
				for (int ct=0; ct<4000; ct++) {
					for (int j=i+1; j<mat.Length; j++) { 
						zeroaverageij(L, mat, i, j); } } 
			i++; } 
			for ( ; i<mat.Length; i+=2) {  
				for (int ct=0; ct<200; ct++) {
					for (int j=i+1; j<mat.Length; j++) { 
						zeroaverageij(L, mat, i, j); } } 
				for (int j=i+2; j<mat.Length; j++) {
					equalizeij(L, mat, i, j); } }           
		} // ===================================================================== 
		void matround(double[][] mat, int n) {
			for (int i=0; i<mat.Length; i++) {
				for (int j=0; j<mat[0].Length; j++) {
					double d = mat[i][j], rd = Math.Round(d, n);
					mat[i][j] = rd; } } 
		} // ===================================================================== 
		int n0 = 2, nlim = 10; 
		double[][][] matn = new double[20][][];
		double[][][] Ln = new double[20][][];
		Font font = new Font("ＭＳ ゴシック", 10); 
		public MyForm() {
			this.ClientSize = new Size(960,980);
			this.Text = "回転の標準化"; 
			for (int n = n0; n<nlim; n++) { // === either even or odd is OK !!! 
				matn[n]=createOrthogonal(n);
				double[][] origcopy=AUtil.newmatcopy(n, n, matn[n]); 
				double[][] ret=AUtil.newmatclear(n, n); 
				double[][] m=AUtil.newmatclear(n, n); 
				Ln[n] = new double[n][];
				for (int i=0; i<n; i++) { Ln[n][i] = new double[n]; 
								Ln[n][i][i] = 1.0; } // === unit matrix !!! 
				getregular(Ln[n], matn[n]); 
				matround(matn[n], 10); }         
			this.Paint += (o,e)=>{ draw(e.Graphics); }; 
			/* 
			Bitmap bm = new Bitmap(this.ClientSize.Width, this.ClientSize.Height);
			draw(bm.getGraphics());  
			for (int i=0; i<bm.Width; i++) { 
				for (int j=0; j<bm.Height; j++) {
					Color c = bm.GetPixel(i, j);
					Color cm = Color.FromArgb(255, c);
					if (c.A <= 10) { cm = Color.FromArgb(0, 0,0,0); } 
					bm.SetPixel(i, j, cm); } } 
			SaveFileDialog dialog = new SaveFileDialog();
			if (dialog.ShowDialog() == DialogResult.OK) {
				string strfile = dialog.FileName;
				bm.Save(strfile); } */ 
		} // =================================================================
		void draw(Graphics g) {
			int J = 0;
			for (int n=n0; n<nlim; n++) { 
				for (int i=0; i<n; i++) {
					for (int j=0; j<n; j++) {
						int x = 20 + j * 100, y = J+20 + i * 18;
						g.DrawString(matn[n][i][j].ToString(),font,
							new SolidBrush(Color.Black),x,y); } } 
				J += 20 + n * 18; }
		} // =================================================================
	} // ##################################################################### 
	public class Ran0 { 
		static ulong IA=16807, IM=2147483647; 
		static protected double AM=1.0/IM; 
		static ulong MASK=123459876; 
		private ulong _idum; 
		protected ulong idum {get {return _idum; } } 
		public Ran0(ulong i) { 
			_idum = i^MASK; 
			_idum %= IM;   
			if (_idum == 0) {_idum=IA; }   } 
		public double step() {  
			kstep(ref _idum, IA, IM); 
			return AM*_idum;   } 
		protected void kstep(ref ulong _idum, ulong IA, ulong IM) { 
			_idum *= IA; 
			_idum %= IM;   } 
	} // #####################################################################