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.

import java.util.*;
import java.text.*;

class Linear {
    private int m, n;
    private double[][] mat;
    private DecimalFormat df;
    
    public Linear(int m,int n) {
        this.m = m; 
        this.n = n;
        mat = new double[m][n+1];
        df = new DecimalFormat("#.00");
    }
    public void getData(Scanner in) {
        for (int i=0;i<m;i++) {
            for (int j=0;j<=n;j++) {
                mat[i][j] = in.nextDouble();
            }
        }
    }
    public void solve() {
        for (int j=0;j<m&&j<n;j++) {
            utilize(j);
            for (int i=j+1;i<m;i++) {
                subtract(i,j);
            }
        }
        int r = rank();
        for (int i=r;i<m;i++) {
            if (!zero(mat[i][n])) {
                System.out.println("NO SOLUTION");
                return;
            }
        }
        for (int i=r-2;i>=0;i--) {
            for (int j=i+1;j<r;j++) {
                subtract(i,j);
            }
        }
        printSol(r);
    }
    private void utilize(int j) {
        int t = m-1;
        for (int i=j;i<=t;i++) {
            if (zero(mat[i][j])) {
                swapRow(i--,t--);
            } else {
                for (int k=n;k>=j;k--) {
                    mat[i][k]/=mat[i][j];
                }
            }
        }
    }
    private void swapRow(int i,int t) {
        for (int j=0;j<=n;j++) {
            double tmp = mat[i][j];
            mat[i][j] = mat[t][j];
            mat[t][j] = tmp;
        }
    }
    private int rank() {
        int r = 0;
        for (;r<m&&r<n;r++) {
            if (zero(mat[r][r])) {
                break;
            }
        }
        return r;
    }
    private void subtract(int i,int j) {
        if (zero(mat[j][j])) {
            return;
        }
        for (int k=n;k>=j;k--) {
            mat[i][k]-=mat[i][j]*mat[j][k];
        }
    }
    private void printSol(int r) {
        System.out.print("X = 	");
        parSol(r);
        for (int j=r;j<n;j++) {
            genSol(r,j);
        }
    }
    private void parSol(int r) {
        System.out.print("(");
        if (r==0) {
            System.out.print(df.format(0));
            for (int i=1;i<n;i++) {
                System.out.print(", "+df.format(0));
            }
        } else {
            System.out.print(df.format(mat[0][n]));
            for (int i=1;i<r;i++) {
                System.out.print(", "+df.format(mat[i][n]));
            }
            for (int i=r;i<n;i++) {
                System.out.print(", "+df.format(0));
            }
        }
        System.out.println(")'");
    }
    private void genSol(int r,int j) {
        System.out.print(" + k"+(j-r+1)+" *	(");
        if (r==0) {
            for (int i=0;i<j;i++) {
                System.out.print(df.format(0)+", ");
            }
            System.out.print(df.format(1));
            for (int i=j+1;i<n;i++) {
                System.out.print(", "+df.format(0));
            }
        } else {
            double len = 1;
            for (int i=0;i<r;i++) {
                len += mat[i][j]*mat[i][j];
            }
            len = Math.sqrt(len);
            System.out.print(df.format(mat[0][j]/len));
            for (int i=1;i<r;i++) {
                System.out.print(", ");
                System.out.print(df.format(mat[i][j]/len));
            }
            for (int i=r;i<j;i++) {
                System.out.print(", "+df.format(0));
            }
            System.out.print(", "+df.format(-1/len));
            for (int i=j+1;i<n;i++) {
                System.out.print(", "+df.format(0));
            }
        }
        System.out.println(")'");
    }
    private boolean zero(double x) {
        return x<0.000001 && x>-0.000001;
    }
}

public class LinearEquation {
    
    public static void main(String[] args) {
        Scanner in = new Scanner(System.in);
        System.out.print("Number of Equations: ");
        int m = in.nextInt();
        System.out.print("Number of Variables: ");
        int n = in.nextInt();
        Linear eq = new Linear(m,n);
        System.out.println("Augmented Matrix:");
        eq.getData(in);
        in.close();
        eq.solve();
    }
}

2. Orthogonal Diagonalization

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

import java.util.*;
import java.text.*;

class Vector {
    private int dim;
    private double[] data;
    
    public Vector(int dim) {
        this.dim = dim;
        data = new double[dim];
    }
    public Vector(Vector other) {
        dim = other.dim;
        data = new double[dim];
        for (int i=0;i<dim;i++) {
            data[i] = other.data[i];
        }
    }
    public int getDim() {
        return dim;
    }
    public double get(int i) {
        if (i<0||i>=dim) {
            throw new RuntimeException("Index out of Bound");
        }
        return data[i];
    }
    public void set(int i,double k) {
        if (i<0||i>=dim) {
            throw new RuntimeException("Index out of Bound");
        }
        data[i] = k;
    }
    public Vector add(Vector other) {
        if (other.dim!=dim) {
            throw new RuntimeException("Invalid Vector Addition");
        }
        Vector val = new Vector(dim);
        for (int i=0;i<dim;i++) {
            val.data[i] = data[i]+other.data[i];
        }
        return val;
    }
    public Vector mult(double k) {
        Vector val = new Vector(dim);
        for (int i=0;i<dim;i++) {
            val.data[i] = data[i]*k;
        }
        return val;
    }
    public double mult(Vector other) {
        if (other.dim!=dim) {
            throw new RuntimeException("Invalid Vector Multiplication");
        }
        double val = 0;
        for (int i=0;i<dim;i++) {
            val += data[i]*other.data[i];
        }
        return val;
    }
    public Vector unify() {
        double len = 0;
        for (int i=0;i<dim;i++) {
            len += data[i]*data[i];
        }
        len = Math.sqrt(len);
        Vector val = new Vector(this);
        if (zero(len)) {
            return val;
        } else if (data[0]<0) {
            len = -len;
        }
        for (int i=0;i<dim;i++) {
            val.data[i] /= len;
        }
        return val;
    }
    private boolean zero(double k) {
        return k>-0.000001&&k<0.000001;
    }
}

class Matrix {
    private int size;
    private double[][] data;
    
    public Matrix(int size) {
        this.size = size;
        data = new double[size][size];
    }
    public Matrix(Matrix other) {
        size = other.size;
        data = new double[size][size];
        for (int i=0;i<size;i++) {
            for (int j=0;j<size;j++) {
                data[i][j] = other.data[i][j];
            }
        }
    }
    public int getSize() {
        return size;
    }
    public double get(int i,int j) {
        if (i<0||i>=size||j<0||j>=size) {
            throw new RuntimeException("Index out of Bound");
        }
        return data[i][j];
    }
    public void set(int i,int j,double k) {
        if (i<0||i>=size||j<0||j>=size) {
            throw new RuntimeException("Index out of Bound");
        }
        data[i][j] = k;
    }
    public double trace() {
        double val = 0;
        for (int i=0;i<size;i++) {
            val += data[i][i];
        }
        return val;
    }
    public double det() {
        double val = 0;
        if (size==1) {
            val += data[0][0];
        } else if (size==2) {
            val += data[0][0]*data[1][1];
            val -= data[0][1]*data[1][0];
        } else if (size==3){
            val += data[0][0]*data[1][1]*data[2][2];
            val += data[0][1]*data[1][2]*data[2][0];
            val += data[1][0]*data[2][1]*data[0][2];
            val -= data[0][2]*data[1][1]*data[2][0];
            val -= data[1][2]*data[2][1]*data[0][0];
            val -= data[2][2]*data[0][1]*data[1][0];
        } else {
            throw new RuntimeException("High-Dim not Supported");
        }
        return val;
    }
    public void getEig(double[] lambda,Vector[] q) {
        // Get Eigenvalues and Eigenvetors:
        if (size==1) {
            lambda[0] = data[0][0];
        } else if (size==2) {
            equSolve(lambda,-trace(),det());
        } else if (size==3) {
            double c = get(0,0)*get(1,1)+get(1,1)*get(2,2)+get(2,2)*get(0,0)
                    - get(0,1)*get(1,0) - get(1,2)*get(2,1) - get(2,0)*get(0,2);
            equSolve(lambda,-trace(),c,-det());
            if (zero(lambda[0]-lambda[2])) {
                double tmp = lambda[2];
                lambda[2] = lambda[1];
                lambda[1] = tmp;
            }
        }
        for (int i=0;i<size;i+=getEigvect(lambda,q,i));
        schmidt(q);
    }
    private int getEigvect(double[] lambda,Vector[] q,int pos) {
        int val = 1;
        for (;pos+val<size;val++) {
            if (!zero(lambda[pos+val]-lambda[pos])) {
                break;
            }
        }
        double[][] mat = new double[size][size];
        for (int i=0;i<size;i++) {
            for (int j=0;j<size;j++) {
                mat[i][j] = -data[i][j];
            }
            mat[i][i] += lambda[pos];
        }
        for (int j=0;j<size;j++) {
            utilize(mat,j);
            for (int i=j+1;i<size;i++) {
                subtract(mat,i,j);
            }
        }
        int r = size-val;
        for (int i=r-2;i>=0;i--) {
            for (int j=i+1;j<r;j++) {
                subtract(mat,i,j);
            }
        }
        for (int j=r;j<size;j++) {
            for (int i=0;i<r;i++) {
                q[pos+j-r].set(i,mat[i][j]);
            }
            q[pos+j-r].set(j,-1);
        }
        return val;
    }
    private void equSolve(double[] rt,double b,double c) {
        // Try to solve a Quadratic Equation: x^2+b*x+c = 0
        double delt = b*b-4*c;
        if (delt<0) {
            throw new RuntimeException("No Real Solution");
        }
        delt = Math.sqrt(delt);
        rt[0] = (-b-delt)/2;
        rt[1] = (-b+delt)/2;
    }
    private void equSolve(double[] rt,double b,double c,double d) {
        // Try to solve a Cubic Equation:  x^3+b*x^2+c*x+d = 0
        double alpha = -b*b*b/27-d/2+b*c/6;
        double beta = c/3-b*b/9;
        if (alpha*alpha+beta*beta*beta>0) {
            throw new RuntimeException("No Real Solution");
        } else if (zero(beta)){
            rt[0] = rt[1] = rt[2] = -b/3; 
            return;
        }
        double tmp = Math.acos(alpha/Math.sqrt(-beta*beta*beta));
        rt[0] = -b/3+2*Math.sqrt(-beta)*Math.cos((tmp-2*Math.PI)/3);
        rt[1] = -b/3+2*Math.sqrt(-beta)*Math.cos(tmp/3);
        rt[2] = -b/3+2*Math.sqrt(-beta)*Math.cos((tmp+2*Math.PI)/3);
    }
    private void schmidt(Vector[] vect) {
        // Gram–Schmidt Orthogonalization
        int num = vect.length;
        Vector[] tmp = new Vector[num];
        for (int i=0;i<num;i++) {
            tmp[i] = new Vector(vect[i]);
        }
        for (int i=0;i<num;i++) {
            for (int j=0;j<i;j++) {
                double k = tmp[i].mult(vect[j]);
                k /= vect[j].mult(vect[j]);
                vect[i] = vect[i].add(vect[j].mult(-k));
            }
            vect[i] = vect[i].unify();
        }
    }
    private void utilize(double[][] mat,int j) {
        int t = size-1;
        for (int i=j;i<=t;i++) {
            if (zero(mat[i][j])) {
                for (int k=j;k<size;k++) {
                    double tmp = mat[i][k];
                    mat[i][k] = mat[t][k];
                    mat[t][k] = tmp;
                }    i --;    t --;
            } else {
                for (int k=size-1;k>=j;k--) {
                    mat[i][k]/=mat[i][j];
                }
            }
        }
    }
    private void subtract(double[][] mat,int i,int j) {
        if (zero(mat[j][j])) {
            return;
        }
        for (int k=size-1;k>=j;k--) {
            mat[i][k] -= mat[i][j]*mat[j][k];
        }
    }
    private boolean zero(double k) {
        return k>-0.000001&&k<0.000001;
    }
}

public class Main {
    public static Scanner in;
    public static DecimalFormat df;
    
    public static void getInput(Matrix coef,Vector rear) {
        int dim = coef.getSize();
        System.out.println("Get the Equation:");
        System.out.print("0 = 	");
        // Second-Degree Terms:
        System.out.print("x^2 * ");
        coef.set(0,0,in.nextDouble());
        System.out.print("  +	");
        System.out.print("y^2 * ");
        coef.set(1,1,in.nextDouble());
        System.out.print("  +	");
        if (dim==3) {
            System.out.print("z^2 * ");
            coef.set(2,2,in.nextDouble());
            System.out.print("  +	");
        }
        System.out.print("x*y * ");
        coef.set(0,1,in.nextDouble()/2);
        coef.set(1,0,coef.get(0,1));
        System.out.print("  +	");
        if (dim==3) {
            System.out.print("y*z * ");
            coef.set(1,2,in.nextDouble()/2);
            coef.set(2,1,coef.get(1,2));
            System.out.print("  +	");
            System.out.print("z*x * ");
            coef.set(2,0,in.nextDouble()/2);
            coef.set(0,2,coef.get(2,0));
            System.out.print("  +	");
        }
        // First-Degree Terms:
        System.out.print("x * ");
        rear.set(0,in.nextDouble());
        System.out.print("  +	");
        System.out.print("y * ");
        rear.set(1,in.nextDouble());
        System.out.print("  +	");
        if (dim==3) {
            System.out.print("z * ");
            rear.set(2,in.nextDouble());
            System.out.print("  +	");
        }
        rear.set(dim,in.nextDouble());
    }
    public static void solve(Matrix coef,Vector rear) {
        int dim = coef.getSize();
        double[] lambda = new double[dim];
        Vector[] q = new Vector[dim];
        for (int i=0;i<dim;i++) {
            q[i] = new Vector(dim);
        }
        coef.getEig(lambda,q);
        printRes(lambda,q,rear);
    }
    public static void printRes(double[] lambda,Vector[] q,Vector rear) {
        int dim = lambda.length;
        double[] bias = new double[dim];
        double cst = rear.get(dim);
        for (int i=0;i<dim;i++) {
            for (int j=0;j<dim;j++) {
                bias[i]+=q[i].get(j)*rear.get(j);
            }
            bias[i] /= (2*lambda[i]);
            cst -= lambda[i]*bias[i]*bias[i];
        }
        System.out.print("0 = 	");
        System.out.print(df.format(lambda[0]));
        System.out.print(" * (");
        translate(0,q,bias[0]);
        System.out.println(")^2");
        System.out.print("  +	");
        System.out.print(df.format(lambda[1]));
        System.out.print(" * (");
        translate(1,q,bias[1]);
        System.out.println(")^2");
        if (dim==3) {
            System.out.print("  +	");
            System.out.print(df.format(lambda[2]));
            System.out.print(" * (");
            translate(2,q,bias[2]);
            System.out.println(")^2");
        }
        System.out.println("  +	"+df.format(cst));
    }
    private static void translate(int pos,Vector[] q,double c) {
        System.out.print(df.format(q[pos].get(0))+"*x");
        System.out.print("+"+df.format(q[pos].get(1))+"*y");
        if (q.length==3) {
            System.out.print("+"+df.format(q[pos].get(2))+"*z");
        }
        if (!zero(c)) {
            System.out.print("+"+df.format(c));
        }
    }
    private static boolean zero(double k) {
        return k>-0.000001&&k<0.000001;
    }
    public static void main(String[] args) {
        in = new Scanner(System.in);
        df = new DecimalFormat("#.0000");
        System.out.println("Curve or Surface?");
        System.out.print("	dim = ");
        int dim = in.nextInt();
        if (dim>=2 && dim<=3) {
            Matrix coef = new Matrix(dim);
            Vector rear = new Vector(dim+1);
            getInput(coef,rear);
            solve(coef,rear);
        } else {
            System.out.println("Sorry, I can't do that.");
        }
        in.close();
    }
}