Linear System and Quadratic Form

　　复习一下线性代数的一些基本概念。

　　一个代数系统指的是定义了一些运算（及其公理）的非空集合。在给定数域内，线性空间是定义在一个非空集合上的代数系统，其上定义了加法和数乘两种运算：(1) 加法满足交换群公理；(2) 数乘满足外结合律并泛化了数域幺元的性质；(3) 两种运算间的外分配律均成立。n维实向量空间、同型实矩阵空间、实函数空间等都是典型的线性空间。
　　给定一个线性空间，我们可以定义线性表出的概念，并由此定义线性相关和线性无关。给定一组线性空间元素，或者说是一个向量组，能线性表出向量组中所有向量的一组线性无关的向量叫做极大线性无关组；极大线性无关组的规模称为原向量组的秩（可证唯一性）。整个线性空间的一个极大线性无关组叫做一组基，线性空间的秩叫做维数。给定一组基，线性空间中的任意向量均可线性表出，并用坐标表示。对于不同的两组基，有一种叫做过渡矩阵的东西可以描述基变换以及同一向量的坐标变换。
　　一个线性空间中，相同运算满足线性空间公理的一个子集称为线性子空间。线性子空间关于交、和这两种运算是封闭的，并且满足一些维数公式。
　　线性映射是定义在两个线性空间之间的同态映射，而矩阵恰恰描述了这种同态关系。矩阵的加法和数乘说明线性映射本身就构成了一个线性空间；矩阵乘法描述了两个同态映射的复合。线性变换是一种特殊的线性映射，描述一个线性空间的自同态。线性变换由方阵描述，其中逆矩阵描述了自同构的逆；不同基下描述同一线性变换的方阵称为相似矩阵（一种等价关系）。
　　给定一组基和一个线性变换，已知象的坐标求原象坐标的过程就是解线性方程组，其通解的结构由系数矩阵和伴随矩阵共同决定。我们定义了一些初等变换以及相应的初等矩阵，左乘初等矩阵等价于初等行变换，右乘初等矩阵等价于初等列变换。相抵是矩阵空间中的一种等价关系，定义为两个矩阵可以由初等变换互相转化，其充要条件是同型且等秩。

　　赋范线性空间是一种特殊的线性空间，其上定义了一种满足正定性、正值齐次性和三角不等式的一元函数，称为范数；Banach空间是按距离完备的赋范线性空间。内积空间是另一种特殊的线性空间，其上定义了一种满足正定性、对称性和双线性性的二元函数，称为内积；由于范数可由内积诱导得到，故内积空间是特殊的赋范线性空间。Hilbert空间是按范数完备的内积空间，也是一种特殊的Banach空间。
　　欧氏空间是一种特殊的实内积空间。在欧氏空间中，给定一个线性无关向量组，我们可以用Gram-Schimdt方法构造一组标准正交基。给定一组基，可以定义一个欧氏空间的度量矩阵；当且仅当基为标准正交基时，度量矩阵为单位阵；不同基下同一欧氏空间的度量矩阵称为合同矩阵（一种等价关系），合同矩阵间满足惯性定律。保持欧氏空间中任意两个向量间内积不变的线性变换叫正交变换，其充要条件是在标准正交基下对应矩阵是正交矩阵。

　　一个n维向量空间上的实值函数对向量的一阶导数称为梯度，二阶导数叫做Hessian矩阵。一个实对称矩阵的特征值与特征向量就是其对应二次型在 $|vec{u}|=1$ 条件下的 Lagrangian 的极值点。一个实对称矩阵的特征值均为实数，特征向量均为实向量，且不同特征值对应的特征向量相互正交；实对称矩阵相似合同于特征值构成的对角矩阵，即 $A=sum_{i=1}^Dlambda_ivec{u}_ivec{u}_i^T$。

1. Gauss-Edmonds Elimination

　　This program can solve a normal linear system and give the general solution if it does exist.

  1 import java.util.*;
  2 import java.text.*;
  3 
  4 class Linear {
  5     private int m, n;
  6     private double[][] mat;
  7     private DecimalFormat df;
  8     
  9     public Linear(int m,int n) {
 10         this.m = m; 
 11         this.n = n;
 12         mat = new double[m][n+1];
 13         df = new DecimalFormat("#.00");
 14     }
 15     public void getData(Scanner in) {
 16         for (int i=0;i<m;i++) {
 17             for (int j=0;j<=n;j++) {
 18                 mat[i][j] = in.nextDouble();
 19             }
 20         }
 21     }
 22     public void solve() {
 23         for (int j=0;j<m&&j<n;j++) {
 24             utilize(j);
 25             for (int i=j+1;i<m;i++) {
 26                 subtract(i,j);
 27             }
 28         }
 29         int r = rank();
 30         for (int i=r;i<m;i++) {
 31             if (!zero(mat[i][n])) {
 32                 System.out.println("NO SOLUTION");
 33                 return;
 34             }
 35         }
 36         for (int i=r-2;i>=0;i--) {
 37             for (int j=i+1;j<r;j++) {
 38                 subtract(i,j);
 39             }
 40         }
 41         printSol(r);
 42     }
 43     private void utilize(int j) {
 44         int t = m-1;
 45         for (int i=j;i<=t;i++) {
 46             if (zero(mat[i][j])) {
 47                 swapRow(i--,t--);
 48             } else {
 49                 for (int k=n;k>=j;k--) {
 50                     mat[i][k]/=mat[i][j];
 51                 }
 52             }
 53         }
 54     }
 55     private void swapRow(int i,int t) {
 56         for (int j=0;j<=n;j++) {
 57             double tmp = mat[i][j];
 58             mat[i][j] = mat[t][j];
 59             mat[t][j] = tmp;
 60         }
 61     }
 62     private int rank() {
 63         int r = 0;
 64         for (;r<m&&r<n;r++) {
 65             if (zero(mat[r][r])) {
 66                 break;
 67             }
 68         }
 69         return r;
 70     }
 71     private void subtract(int i,int j) {
 72         if (zero(mat[j][j])) {
 73             return;
 74         }
 75         for (int k=n;k>=j;k--) {
 76             mat[i][k]-=mat[i][j]*mat[j][k];
 77         }
 78     }
 79     private void printSol(int r) {
 80         System.out.print("X = 	");
 81         parSol(r);
 82         for (int j=r;j<n;j++) {
 83             genSol(r,j);
 84         }
 85     }
 86     private void parSol(int r) {
 87         System.out.print("(");
 88         if (r==0) {
 89             System.out.print(df.format(0));
 90             for (int i=1;i<n;i++) {
 91                 System.out.print(", "+df.format(0));
 92             }
 93         } else {
 94             System.out.print(df.format(mat[0][n]));
 95             for (int i=1;i<r;i++) {
 96                 System.out.print(", "+df.format(mat[i][n]));
 97             }
 98             for (int i=r;i<n;i++) {
 99                 System.out.print(", "+df.format(0));
100             }
101         }
102         System.out.println(")'");
103     }
104     private void genSol(int r,int j) {
105         System.out.print(" + k"+(j-r+1)+" *	(");
106         if (r==0) {
107             for (int i=0;i<j;i++) {
108                 System.out.print(df.format(0)+", ");
109             }
110             System.out.print(df.format(1));
111             for (int i=j+1;i<n;i++) {
112                 System.out.print(", "+df.format(0));
113             }
114         } else {
115             double len = 1;
116             for (int i=0;i<r;i++) {
117                 len += mat[i][j]*mat[i][j];
118             }
119             len = Math.sqrt(len);
120             System.out.print(df.format(mat[0][j]/len));
121             for (int i=1;i<r;i++) {
122                 System.out.print(", ");
123                 System.out.print(df.format(mat[i][j]/len));
124             }
125             for (int i=r;i<j;i++) {
126                 System.out.print(", "+df.format(0));
127             }
128             System.out.print(", "+df.format(-1/len));
129             for (int i=j+1;i<n;i++) {
130                 System.out.print(", "+df.format(0));
131             }
132         }
133         System.out.println(")'");
134     }
135     private boolean zero(double x) {
136         return x<0.000001 && x>-0.000001;
137     }
138 }
139 
140 public class LinearEquation {
141     
142     public static void main(String[] args) {
143         Scanner in = new Scanner(System.in);
144         System.out.print("Number of Equations: ");
145         int m = in.nextInt();
146         System.out.print("Number of Variables: ");
147         int n = in.nextInt();
148         Linear eq = new Linear(m,n);
149         System.out.println("Augmented Matrix:");
150         eq.getData(in);
151         in.close();
152         eq.solve();
153     }
154 }

2. Orthogonal Diagonalization

　　This program can show the standard equation of a 2-dimensional quadratic curve or a 3-dimensional quadratic surface:

  1 import java.util.*;
  2 import java.text.*;
  3 
  4 class Vector {
  5     private int dim;
  6     private double[] data;
  7     
  8     public Vector(int dim) {
  9         this.dim = dim;
 10         data = new double[dim];
 11     }
 12     public Vector(Vector other) {
 13         dim = other.dim;
 14         data = new double[dim];
 15         for (int i=0;i<dim;i++) {
 16             data[i] = other.data[i];
 17         }
 18     }
 19     public int getDim() {
 20         return dim;
 21     }
 22     public double get(int i) {
 23         if (i<0||i>=dim) {
 24             throw new RuntimeException("Index out of Bound");
 25         }
 26         return data[i];
 27     }
 28     public void set(int i,double k) {
 29         if (i<0||i>=dim) {
 30             throw new RuntimeException("Index out of Bound");
 31         }
 32         data[i] = k;
 33     }
 34     public Vector add(Vector other) {
 35         if (other.dim!=dim) {
 36             throw new RuntimeException("Invalid Vector Addition");
 37         }
 38         Vector val = new Vector(dim);
 39         for (int i=0;i<dim;i++) {
 40             val.data[i] = data[i]+other.data[i];
 41         }
 42         return val;
 43     }
 44     public Vector mult(double k) {
 45         Vector val = new Vector(dim);
 46         for (int i=0;i<dim;i++) {
 47             val.data[i] = data[i]*k;
 48         }
 49         return val;
 50     }
 51     public double mult(Vector other) {
 52         if (other.dim!=dim) {
 53             throw new RuntimeException("Invalid Vector Multiplication");
 54         }
 55         double val = 0;
 56         for (int i=0;i<dim;i++) {
 57             val += data[i]*other.data[i];
 58         }
 59         return val;
 60     }
 61     public Vector unify() {
 62         double len = 0;
 63         for (int i=0;i<dim;i++) {
 64             len += data[i]*data[i];
 65         }
 66         len = Math.sqrt(len);
 67         Vector val = new Vector(this);
 68         if (zero(len)) {
 69             return val;
 70         } else if (data[0]<0) {
 71             len = -len;
 72         }
 73         for (int i=0;i<dim;i++) {
 74             val.data[i] /= len;
 75         }
 76         return val;
 77     }
 78     private boolean zero(double k) {
 79         return k>-0.000001&&k<0.000001;
 80     }
 81 }
 82 
 83 class Matrix {
 84     private int size;
 85     private double[][] data;
 86     
 87     public Matrix(int size) {
 88         this.size = size;
 89         data = new double[size][size];
 90     }
 91     public Matrix(Matrix other) {
 92         size = other.size;
 93         data = new double[size][size];
 94         for (int i=0;i<size;i++) {
 95             for (int j=0;j<size;j++) {
 96                 data[i][j] = other.data[i][j];
 97             }
 98         }
 99     }
100     public int getSize() {
101         return size;
102     }
103     public double get(int i,int j) {
104         if (i<0||i>=size||j<0||j>=size) {
105             throw new RuntimeException("Index out of Bound");
106         }
107         return data[i][j];
108     }
109     public void set(int i,int j,double k) {
110         if (i<0||i>=size||j<0||j>=size) {
111             throw new RuntimeException("Index out of Bound");
112         }
113         data[i][j] = k;
114     }
115     public double trace() {
116         double val = 0;
117         for (int i=0;i<size;i++) {
118             val += data[i][i];
119         }
120         return val;
121     }
122     public double det() {
123         double val = 0;
124         if (size==1) {
125             val += data[0][0];
126         } else if (size==2) {
127             val += data[0][0]*data[1][1];
128             val -= data[0][1]*data[1][0];
129         } else if (size==3){
130             val += data[0][0]*data[1][1]*data[2][2];
131             val += data[0][1]*data[1][2]*data[2][0];
132             val += data[1][0]*data[2][1]*data[0][2];
133             val -= data[0][2]*data[1][1]*data[2][0];
134             val -= data[1][2]*data[2][1]*data[0][0];
135             val -= data[2][2]*data[0][1]*data[1][0];
136         } else {
137             throw new RuntimeException("High-Dim not Supported");
138         }
139         return val;
140     }
141     public void getEig(double[] lambda,Vector[] q) {
142         // Get Eigenvalues and Eigenvetors:
143         if (size==1) {
144             lambda[0] = data[0][0];
145         } else if (size==2) {
146             equSolve(lambda,-trace(),det());
147         } else if (size==3) {
148             double c = get(0,0)*get(1,1)+get(1,1)*get(2,2)+get(2,2)*get(0,0)
149                     - get(0,1)*get(1,0) - get(1,2)*get(2,1) - get(2,0)*get(0,2);
150             equSolve(lambda,-trace(),c,-det());
151             if (zero(lambda[0]-lambda[2])) {
152                 double tmp = lambda[2];
153                 lambda[2] = lambda[1];
154                 lambda[1] = tmp;
155             }
156         }
157         for (int i=0;i<size;i+=getEigvect(lambda,q,i));
158         schmidt(q);
159     }
160     private int getEigvect(double[] lambda,Vector[] q,int pos) {
161         int val = 1;
162         for (;pos+val<size;val++) {
163             if (!zero(lambda[pos+val]-lambda[pos])) {
164                 break;
165             }
166         }
167         double[][] mat = new double[size][size];
168         for (int i=0;i<size;i++) {
169             for (int j=0;j<size;j++) {
170                 mat[i][j] = -data[i][j];
171             }
172             mat[i][i] += lambda[pos];
173         }
174         for (int j=0;j<size;j++) {
175             utilize(mat,j);
176             for (int i=j+1;i<size;i++) {
177                 subtract(mat,i,j);
178             }
179         }
180         int r = size-val;
181         for (int i=r-2;i>=0;i--) {
182             for (int j=i+1;j<r;j++) {
183                 subtract(mat,i,j);
184             }
185         }
186         for (int j=r;j<size;j++) {
187             for (int i=0;i<r;i++) {
188                 q[pos+j-r].set(i,mat[i][j]);
189             }
190             q[pos+j-r].set(j,-1);
191         }
192         return val;
193     }
194     private void equSolve(double[] rt,double b,double c) {
195         // Try to solve a Quadratic Equation: x^2+b*x+c = 0
196         double delt = b*b-4*c;
197         if (delt<0) {
198             throw new RuntimeException("No Real Solution");
199         }
200         delt = Math.sqrt(delt);
201         rt[0] = (-b-delt)/2;
202         rt[1] = (-b+delt)/2;
203     }
204     private void equSolve(double[] rt,double b,double c,double d) {
205         // Try to solve a Cubic Equation:  x^3+b*x^2+c*x+d = 0
206         double alpha = -b*b*b/27-d/2+b*c/6;
207         double beta = c/3-b*b/9;
208         if (alpha*alpha+beta*beta*beta>0) {
209             throw new RuntimeException("No Real Solution");
210         } else if (zero(beta)){
211             rt[0] = rt[1] = rt[2] = -b/3; 
212             return;
213         }
214         double tmp = Math.acos(alpha/Math.sqrt(-beta*beta*beta));
215         rt[0] = -b/3+2*Math.sqrt(-beta)*Math.cos((tmp-2*Math.PI)/3);
216         rt[1] = -b/3+2*Math.sqrt(-beta)*Math.cos(tmp/3);
217         rt[2] = -b/3+2*Math.sqrt(-beta)*Math.cos((tmp+2*Math.PI)/3);
218     }
219     private void schmidt(Vector[] vect) {
220         // Gram–Schmidt Orthogonalization
221         int num = vect.length;
222         Vector[] tmp = new Vector[num];
223         for (int i=0;i<num;i++) {
224             tmp[i] = new Vector(vect[i]);
225         }
226         for (int i=0;i<num;i++) {
227             for (int j=0;j<i;j++) {
228                 double k = tmp[i].mult(vect[j]);
229                 k /= vect[j].mult(vect[j]);
230                 vect[i] = vect[i].add(vect[j].mult(-k));
231             }
232             vect[i] = vect[i].unify();
233         }
234     }
235     private void utilize(double[][] mat,int j) {
236         int t = size-1;
237         for (int i=j;i<=t;i++) {
238             if (zero(mat[i][j])) {
239                 for (int k=j;k<size;k++) {
240                     double tmp = mat[i][k];
241                     mat[i][k] = mat[t][k];
242                     mat[t][k] = tmp;
243                 }    i --;    t --;
244             } else {
245                 for (int k=size-1;k>=j;k--) {
246                     mat[i][k]/=mat[i][j];
247                 }
248             }
249         }
250     }
251     private void subtract(double[][] mat,int i,int j) {
252         if (zero(mat[j][j])) {
253             return;
254         }
255         for (int k=size-1;k>=j;k--) {
256             mat[i][k] -= mat[i][j]*mat[j][k];
257         }
258     }
259     private boolean zero(double k) {
260         return k>-0.000001&&k<0.000001;
261     }
262 }
263 
264 public class Main {
265     public static Scanner in;
266     public static DecimalFormat df;
267     
268     public static void getInput(Matrix coef,Vector rear) {
269         int dim = coef.getSize();
270         System.out.println("Get the Equation:");
271         System.out.print("0 = 	");
272         // Second-Degree Terms:
273         System.out.print("x^2 * ");
274         coef.set(0,0,in.nextDouble());
275         System.out.print("  +	");
276         System.out.print("y^2 * ");
277         coef.set(1,1,in.nextDouble());
278         System.out.print("  +	");
279         if (dim==3) {
280             System.out.print("z^2 * ");
281             coef.set(2,2,in.nextDouble());
282             System.out.print("  +	");
283         }
284         System.out.print("x*y * ");
285         coef.set(0,1,in.nextDouble()/2);
286         coef.set(1,0,coef.get(0,1));
287         System.out.print("  +	");
288         if (dim==3) {
289             System.out.print("y*z * ");
290             coef.set(1,2,in.nextDouble()/2);
291             coef.set(2,1,coef.get(1,2));
292             System.out.print("  +	");
293             System.out.print("z*x * ");
294             coef.set(2,0,in.nextDouble()/2);
295             coef.set(0,2,coef.get(2,0));
296             System.out.print("  +	");
297         }
298         // First-Degree Terms:
299         System.out.print("x * ");
300         rear.set(0,in.nextDouble());
301         System.out.print("  +	");
302         System.out.print("y * ");
303         rear.set(1,in.nextDouble());
304         System.out.print("  +	");
305         if (dim==3) {
306             System.out.print("z * ");
307             rear.set(2,in.nextDouble());
308             System.out.print("  +	");
309         }
310         rear.set(dim,in.nextDouble());
311     }
312     public static void solve(Matrix coef,Vector rear) {
313         int dim = coef.getSize();
314         double[] lambda = new double[dim];
315         Vector[] q = new Vector[dim];
316         for (int i=0;i<dim;i++) {
317             q[i] = new Vector(dim);
318         }
319         coef.getEig(lambda,q);
320         printRes(lambda,q,rear);
321     }
322     public static void printRes(double[] lambda,Vector[] q,Vector rear) {
323         int dim = lambda.length;
324         double[] bias = new double[dim];
325         double cst = rear.get(dim);
326         for (int i=0;i<dim;i++) {
327             for (int j=0;j<dim;j++) {
328                 bias[i]+=q[i].get(j)*rear.get(j);
329             }
330             bias[i] /= (2*lambda[i]);
331             cst -= lambda[i]*bias[i]*bias[i];
332         }
333         System.out.print("0 = 	");
334         System.out.print(df.format(lambda[0]));
335         System.out.print(" * (");
336         translate(0,q,bias[0]);
337         System.out.println(")^2");
338         System.out.print("  +	");
339         System.out.print(df.format(lambda[1]));
340         System.out.print(" * (");
341         translate(1,q,bias[1]);
342         System.out.println(")^2");
343         if (dim==3) {
344             System.out.print("  +	");
345             System.out.print(df.format(lambda[2]));
346             System.out.print(" * (");
347             translate(2,q,bias[2]);
348             System.out.println(")^2");
349         }
350         System.out.println("  +	"+df.format(cst));
351     }
352     private static void translate(int pos,Vector[] q,double c) {
353         System.out.print(df.format(q[pos].get(0))+"*x");
354         System.out.print("+"+df.format(q[pos].get(1))+"*y");
355         if (q.length==3) {
356             System.out.print("+"+df.format(q[pos].get(2))+"*z");
357         }
358         if (!zero(c)) {
359             System.out.print("+"+df.format(c));
360         }
361     }
362     private static boolean zero(double k) {
363         return k>-0.000001&&k<0.000001;
364     }
365     public static void main(String[] args) {
366         in = new Scanner(System.in);
367         df = new DecimalFormat("#.0000");
368         System.out.println("Curve or Surface?");
369         System.out.print("	dim = ");
370         int dim = in.nextInt();
371         if (dim>=2 && dim<=3) {
372             Matrix coef = new Matrix(dim);
373             Vector rear = new Vector(dim+1);
374             getInput(coef,rear);
375             solve(coef,rear);
376         } else {
377             System.out.println("Sorry, I can't do that.");
378         }
379         in.close();
380     }
381 }