QR decomposition: Difference between revisions
Content added Content deleted
m (→{{header|ATS}}) |
(New post which completes the task without using external libraries. In addition to 4 existing posts, none of which complete the task, and use the external libraries "JAMA", "Colt", "Apache Commons Math" and "la4j".) |
||
Line 2,950: | Line 2,950: | ||
0,000 -175,000 70,000 |
0,000 -175,000 70,000 |
||
0,000 0,000 35,000</pre> |
0,000 0,000 35,000</pre> |
||
===Without external libraries=== |
|||
<syntaxhighlight lang="java"> |
|||
import java.util.ArrayList; |
|||
import java.util.List; |
|||
public final class QRDecomposition { |
|||
public static void main(String[] aArgs) { |
|||
final double[][] data = new double [][] { { 12.0, -51.0, 4.0 }, |
|||
{ 6.0, 167.0, -68.0 }, |
|||
{ -4.0, 24.0, -41.0 }, |
|||
{ -1.0, 1.0, 0.0 }, |
|||
{ 2.0, 0.0, 3.0 } }; |
|||
// Task 1 |
|||
Matrix A = new Matrix(data); |
|||
A.display("Initial matrix A:"); |
|||
MatrixPair pair = householder(A); |
|||
Matrix Q = pair.q; |
|||
Matrix R = pair.r; |
|||
Q.display("Matrix Q:"); |
|||
R.display("Matrix R:"); |
|||
Matrix result = Q.multiply(R); |
|||
result.display("Matrix Q * R:"); |
|||
// Task 2 |
|||
Matrix x = new Matrix ( new double[][] { { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0 } } ); |
|||
Matrix y = new Matrix( |
|||
new double[][] { { 1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0 } } ); |
|||
result = fitPolynomial(x, y, 2); |
|||
result.display("Result of fitting polynomial:"); |
|||
} |
|||
private static MatrixPair householder(Matrix aMatrix) { |
|||
final int rowCount = aMatrix.getRowCount(); |
|||
final int columnCount = aMatrix.getColumnCount(); |
|||
List<Matrix> versionsOfQ = new ArrayList<Matrix>(rowCount); |
|||
Matrix z = new Matrix(aMatrix); |
|||
Matrix z1 = new Matrix(rowCount, columnCount); |
|||
for ( int k = 0; k < columnCount && k < rowCount - 1; k++ ) { |
|||
Matrix vectorE = new Matrix( new double[rowCount][1] ); |
|||
z1 = z.minor(k); |
|||
Matrix vectorX = z1.column(k); |
|||
double magnitudeX = vectorX.magnitude(); |
|||
if ( aMatrix.getEntry(k, k) > 0 ) { |
|||
magnitudeX = -magnitudeX; |
|||
} |
|||
for ( int i = 0; i < vectorE.size(); i++ ) { |
|||
vectorE.setEntry(i, 0, ( i == k ) ? 1 : 0); |
|||
} |
|||
vectorE = vectorE.scalarMultiply(magnitudeX).add(vectorX); |
|||
vectorE = vectorE.unit(); |
|||
versionsOfQ.add(householderFactor(vectorE)); |
|||
z = versionsOfQ.get(k).multiply(z1); |
|||
} |
|||
Matrix Q = versionsOfQ.get(0); |
|||
for ( int i = 1; i < columnCount && i < rowCount - 1; i++ ) { |
|||
Q = versionsOfQ.get(i).multiply(Q); |
|||
} |
|||
Matrix R = Q.multiply(aMatrix); |
|||
Q = Q.transpose(); |
|||
return new MatrixPair(R, Q); |
|||
} |
|||
public static Matrix householderFactor(Matrix aVector) { |
|||
if ( aVector.getColumnCount() != 1 ) { |
|||
throw new RuntimeException("Incompatible matrix dimensions."); |
|||
} |
|||
final int size = aVector.size(); |
|||
double[][] newData = new double[size][size]; |
|||
for ( int i = 0; i < size; i++ ) { |
|||
for ( int j = 0; j < size; j++ ) { |
|||
newData[i][j] = -2 * aVector.getEntry(i, 0) * aVector.getEntry(j, 0); |
|||
} |
|||
} |
|||
for ( int i = 0; i < size; i++ ) { |
|||
newData[i][i] += 1; |
|||
} |
|||
return new Matrix(newData); |
|||
} |
|||
private static Matrix fitPolynomial(Matrix aX, Matrix aY, int aPolynomialDegree) { |
|||
Matrix vandermonde = new Matrix(aX.getColumnCount(), aPolynomialDegree + 1); |
|||
for ( int i = 0; i < aX.getColumnCount(); i++ ) { |
|||
for ( int j = 0; j < aPolynomialDegree + 1; j++ ) { |
|||
vandermonde.setEntry(i, j, Math.pow(aX.getEntry(0, i), j)); |
|||
} |
|||
} |
|||
return leastSquares(vandermonde, aY.transpose()); |
|||
} |
|||
private static Matrix leastSquares(Matrix aVandermonde, Matrix aB) { |
|||
MatrixPair pair = householder(aVandermonde); |
|||
return solveUpperTriangular(pair.r, pair.q.transpose().multiply(aB)); |
|||
} |
|||
private static Matrix solveUpperTriangular(Matrix aR, Matrix aB) { |
|||
final int columnCount = aR.getColumnCount(); |
|||
Matrix result = new Matrix(columnCount, 1); |
|||
for ( int k = columnCount - 1; k >= 0; k-- ) { |
|||
double total = 0.0; |
|||
for ( int j = k + 1; j < columnCount; j++ ) { |
|||
total += aR.getEntry(k, j) * result.getEntry(j, 0); |
|||
} |
|||
result.setEntry(k, 0, ( aB.getEntry(k, 0) - total ) / aR.getEntry(k, k)); |
|||
} |
|||
return result; |
|||
} |
|||
private static record MatrixPair(Matrix r, Matrix q) {} |
|||
} |
|||
final class Matrix { |
|||
public Matrix(double[][] aData) { |
|||
rowCount = aData.length; |
|||
columnCount = aData[0].length; |
|||
data = new double[rowCount][columnCount]; |
|||
for ( int i = 0; i < rowCount; i++ ) { |
|||
for ( int j = 0; j < columnCount; j++ ) { |
|||
data[i][j] = aData[i][j]; |
|||
} |
|||
} |
|||
} |
|||
public Matrix(Matrix aMatrix) { |
|||
this(aMatrix.data); |
|||
} |
|||
public Matrix(int aRows, int aCols) { |
|||
this( new double[aRows][aCols] ); |
|||
} |
|||
public Matrix add(Matrix aOther) { |
|||
if ( aOther.rowCount != rowCount || aOther.columnCount != columnCount ) { |
|||
throw new IllegalArgumentException("Incompatible matrix dimensions."); |
|||
} |
|||
Matrix result = new Matrix(data); |
|||
for ( int i = 0; i < rowCount; i++ ) { |
|||
for ( int j = 0; j < columnCount; j++ ) { |
|||
result.data[i][j] = data[i][j] + aOther.data[i][j]; |
|||
} |
|||
} |
|||
return result; |
|||
} |
|||
public Matrix multiply(Matrix aOther) { |
|||
if ( columnCount != aOther.rowCount ) { |
|||
throw new IllegalArgumentException("Incompatible matrix dimensions."); |
|||
} |
|||
Matrix result = new Matrix(rowCount, aOther.columnCount); |
|||
for ( int i = 0; i < rowCount; i++ ) { |
|||
for ( int j = 0; j < aOther.columnCount; j++ ) { |
|||
for ( int k = 0; k < rowCount; k++ ) { |
|||
result.data[i][j] += data[i][k] * aOther.data[k][j]; |
|||
} |
|||
} |
|||
} |
|||
return result; |
|||
} |
|||
public Matrix transpose() { |
|||
Matrix result = new Matrix(columnCount, rowCount); |
|||
for ( int i = 0; i < rowCount; i++ ) { |
|||
for ( int j = 0; j < columnCount; j++ ) { |
|||
result.data[j][i] = data[i][j]; |
|||
} |
|||
} |
|||
return result; |
|||
} |
|||
public Matrix minor(int aIndex) { |
|||
Matrix result = new Matrix(rowCount, columnCount); |
|||
for ( int i = 0; i < aIndex; i++ ) { |
|||
result.setEntry(i, i, 1.0); |
|||
} |
|||
for ( int i = aIndex; i < rowCount; i++ ) { |
|||
for ( int j = aIndex; j < columnCount; j++ ) { |
|||
result.setEntry(i, j, data[i][j]); |
|||
} |
|||
} |
|||
return result; |
|||
} |
|||
public Matrix column(int aColumn) { |
|||
Matrix result = new Matrix(rowCount, 1); |
|||
for ( int i = 0; i < rowCount; i++ ) { |
|||
result.setEntry(i, 0, data[i][aColumn]); |
|||
} |
|||
return result; |
|||
} |
|||
public Matrix scalarMultiply(double aValue) { |
|||
if ( columnCount != 1 ) { |
|||
throw new IllegalArgumentException("Incompatible matrix dimension."); |
|||
} |
|||
Matrix result = new Matrix(rowCount, columnCount); |
|||
for ( int i = 0; i < data.length; i++ ) { |
|||
result.data[i][0] = data[i][0] * aValue; |
|||
} |
|||
return result; |
|||
} |
|||
public Matrix unit() { |
|||
if ( columnCount != 1 ) { |
|||
throw new IllegalArgumentException("Incompatible matrix dimensions."); |
|||
} |
|||
final double magnitude = magnitude(); |
|||
Matrix result = new Matrix(rowCount, columnCount); |
|||
for ( int i = 0; i < data.length; i++ ) { |
|||
result.data[i][0] = data[i][0] / magnitude; |
|||
} |
|||
return result; |
|||
} |
|||
public double magnitude() { |
|||
if ( columnCount != 1 ) { |
|||
throw new IllegalArgumentException("Incompatible matrix dimensions."); |
|||
} |
|||
double norm = 0.0; |
|||
for ( int i = 0; i < data.length; i++ ) { |
|||
norm += data[i][0] * data[i][0]; |
|||
} |
|||
return Math.sqrt(norm); |
|||
} |
|||
public int size() { |
|||
if ( columnCount != 1 ) { |
|||
throw new IllegalArgumentException("Incompatible matrix dimensions."); |
|||
} |
|||
return data.length; |
|||
} |
|||
public void display(String aTitle) { |
|||
System.out.println(aTitle); |
|||
for ( int i = 0; i < rowCount; i++ ) { |
|||
for ( int j = 0; j < columnCount; j++ ) { |
|||
System.out.print(String.format("%9.4f", data[i][j])); |
|||
} |
|||
System.out.println(); |
|||
} |
|||
System.out.println(); |
|||
} |
|||
public double getEntry(int aRow, int aCol) { |
|||
return data[aRow][aCol]; |
|||
} |
|||
public void setEntry(int aRow, int aCol, double aValue) { |
|||
data[aRow][aCol] = aValue; |
|||
} |
|||
public int getRowCount() { |
|||
return rowCount; |
|||
} |
|||
public int getColumnCount() { |
|||
return columnCount; |
|||
} |
|||
private final int rowCount; |
|||
private final int columnCount; |
|||
private final double[][] data; |
|||
} |
|||
</syntaxhighlight> |
|||
{{ out }} |
|||
<pre> |
|||
Initial matrix A: |
|||
12.0000 -51.0000 4.0000 |
|||
6.0000 167.0000 -68.0000 |
|||
-4.0000 24.0000 -41.0000 |
|||
-1.0000 1.0000 0.0000 |
|||
2.0000 0.0000 3.0000 |
|||
Matrix Q: |
|||
0.8464 -0.3913 0.3431 0.0815 0.0781 |
|||
0.4232 0.9041 -0.0293 0.0258 0.0447 |
|||
-0.2821 0.1704 0.9329 -0.0474 -0.1374 |
|||
-0.0705 0.0140 -0.0011 0.9804 -0.1836 |
|||
0.1411 -0.0167 -0.1058 -0.1713 -0.9692 |
|||
Matrix R: |
|||
14.1774 20.6666 -13.4016 |
|||
-0.0000 175.0425 -70.0803 |
|||
0.0000 0.0000 -35.2015 |
|||
-0.0000 -0.0000 -0.0000 |
|||
0.0000 0.0000 -0.0000 |
|||
Matrix Q * R: |
|||
12.0000 -51.0000 4.0000 |
|||
6.0000 167.0000 -68.0000 |
|||
-4.0000 24.0000 -41.0000 |
|||
-1.0000 1.0000 -0.0000 |
|||
2.0000 -0.0000 3.0000 |
|||
Result of fitting polynomial: |
|||
1.0000 |
|||
2.0000 |
|||
3.0000 |
|||
</pre> |
|||
=={{header|jq}}== |
=={{header|jq}}== |