Linear System and Quadratic Form
复习一下线性代数的一些基本概念。
一个代数系统指的是定义了一些运算（及其公理）的非空集合。在给定数域内，线性空间是定义在一个非空集合上的代数系统，其上定义了加法和数乘两种运算：(1) 加法满足交换群公理；(2) 数乘满足外结合律并泛化了数域幺元的性质；(3) 两种运算间的外分配律均成立。n维实向量空间、同型实矩阵空间、实函数空间等都是典型的线性空间。
给定一个线性空间，我们可以定义线性表出的概念，并由此定义线性相关和线性无关。给定一组线性空间元素，或者说是一个向量组，能线性表出向量组中所有向量的一组线性无关的向量叫做极大线性无关组；极大线性无关组的规模称为原向量组的秩（可证唯一性）。整个线性空间的一个极大线性无关组叫做一组基，线性空间的秩叫做维数。给定一组基，线性空间中的任意向量均可线性表出，并用坐标表示。对于不同的两组基，有一种叫做过渡矩阵的东西可以描述基变换以及同一向量的坐标变换。
一个线性空间中，相同运算满足线性空间公理的一个子集称为线性子空间。线性子空间关于交、和这两种运算是封闭的，并且满足一些维数公式。
线性映射是定义在两个线性空间之间的同态映射，而矩阵恰恰描述了这种同态关系。矩阵的加法和数乘说明线性映射本身就构成了一个线性空间；矩阵乘法描述了两个同态映射的复合。线性变换是一种特殊的线性映射，描述一个线性空间的自同态。线性变换由方阵描述，其中逆矩阵描述了自同构的逆；不同基下描述同一线性变换的方阵称为相似矩阵（一种等价关系）。
给定一组基和一个线性变换，已知象的坐标求原象坐标的过程就是解线性方程组，其通解的结构由系数矩阵和伴随矩阵共同决定。我们定义了一些初等变换以及相应的初等矩阵，左乘初等矩阵等价于初等行变换，右乘初等矩阵等价于初等列变换。相抵是矩阵空间中的一种等价关系，定义为两个矩阵可以由初等变换互相转化，其充要条件是同型且等秩。
赋范线性空间是一种特殊的线性空间，其上定义了一种满足正定性、正值齐次性和三角不等式的一元函数，称为范数；Banach空间是按距离完备的赋范线性空间。内积空间是另一种特殊的线性空间，其上定义了一种满足正定性、对称性和双线性性的二元函数，称为内积；由于范数可由内积诱导得到，故内积空间是特殊的赋范线性空间。Hilbert空间是按范数完备的内积空间，也是一种特殊的Banach空间。
欧氏空间是一种特殊的实内积空间。在欧氏空间中，给定一个线性无关向量组，我们可以用Gram-Schimdt方法构造一组标准正交基。给定一组基，可以定义一个欧氏空间的度量矩阵；当且仅当基为标准正交基时，度量矩阵为单位阵；不同基下同一欧氏空间的度量矩阵称为合同矩阵（一种等价关系），合同矩阵间满足惯性定律。保持欧氏空间中任意两个向量间内积不变的线性变换叫正交变换，其充要条件是在标准正交基下对应矩阵是正交矩阵。
一个n维向量空间上的实值函数对向量的一阶导数称为梯度，二阶导数叫做Hessian矩阵。一个实对称矩阵的特征值与特征向量就是其对应二次型在 $|\vec{u}|=1$ 条件下的 Lagrangian 的极值点。一个实对称矩阵的特征值均为实数，特征向量均为实向量，且不同特征值对应的特征向量相互正交；实对称矩阵相似合同于特征值构成的对角矩阵，即 $A=\sum_{i=1}^D\lambda_i\vec{u}_i\vec{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 = \t"); 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)+" *\t("); 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 = \t"); 272 // Second-Degree Terms: 273 System.out.print("x^2 * "); 274 coef.set(0,0,in.nextDouble()); 275 System.out.print(" +\t"); 276 System.out.print("y^2 * "); 277 coef.set(1,1,in.nextDouble()); 278 System.out.print(" +\t"); 279 if (dim==3) { 280 System.out.print("z^2 * "); 281 coef.set(2,2,in.nextDouble()); 282 System.out.print(" +\t"); 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(" +\t"); 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(" +\t"); 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(" +\t"); 297 } 298 // First-Degree Terms: 299 System.out.print("x * "); 300 rear.set(0,in.nextDouble()); 301 System.out.print(" +\t"); 302 System.out.print("y * "); 303 rear.set(1,in.nextDouble()); 304 System.out.print(" +\t"); 305 if (dim==3) { 306 System.out.print("z * "); 307 rear.set(2,in.nextDouble()); 308 System.out.print(" +\t"); 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 = \t"); 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(" +\t"); 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(" +\t"); 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(" +\t"+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("\tdim = "); 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 }
