Rosetta Code

QR decomposition: Difference between revisions

Page Discussion

QR decomposition (view source)

Revision as of 12:03, 29 January 2024

17,409 bytes added ,  4 months ago
m
→‎{{header|Wren}}: Minor tidy
m (→‎{{header|Wren}}: Minor tidy)
 
(5 intermediate revisions by one other user not shown)
Line 1,951:
-1.000 1.000 -0.000
2.000 -0.000 3.000
</pre>
 
===With Polynomial Fitting===
<syntaxhighlight lang="c++">
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>
 
class Matrix {
public:
Matrix(const std::vector<std::vector<double>>& data) : data(data) {
initialise();
}
 
Matrix(const Matrix& matrix) : data(matrix.data) {
initialise();
}
 
Matrix(const uint64_t& row_count, const uint64_t& column_count) {
data.assign(row_count, std::vector<double>(column_count, 0.0));
initialise();
}
 
Matrix add(const Matrix& other) {
if ( other.row_count != row_count || other.column_count != column_count ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
Matrix result(data);
for ( int32_t i = 0; i < row_count; ++i ) {
for ( int32_t j = 0; j < column_count; ++j ) {
result.data[i][j] = data[i][j] + other.data[i][j];
}
}
return result;
}
 
Matrix multiply(const Matrix& other) {
if ( column_count != other.row_count ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
Matrix result(row_count, other.column_count);
for ( int32_t i = 0; i < row_count; ++i ) {
for ( int32_t j = 0; j < other.column_count; ++j ) {
for ( int32_t k = 0; k < row_count; k++ ) {
result.data[i][j] += data[i][k] * other.data[k][j];
}
}
}
return result;
}
 
Matrix transpose() {
Matrix result(column_count, row_count);
for ( int32_t i = 0; i < row_count; ++i ) {
for ( int32_t j = 0; j < column_count; ++j ) {
result.data[j][i] = data[i][j];
}
}
return result;
}
 
Matrix minor(const int32_t& index) {
Matrix result(row_count, column_count);
for ( int32_t i = 0; i < index; ++i ) {
result.set_entry(i, i, 1.0);
}
 
for ( int32_t i = index; i < row_count; ++i ) {
for ( int32_t j = index; j < column_count; ++j ) {
result.set_entry(i, j, data[i][j]);
}
}
return result;
}
 
Matrix column(const int32_t& index) {
Matrix result(row_count, 1);
for ( int32_t i = 0; i < row_count; ++i ) {
result.set_entry(i, 0, data[i][index]);
}
return result;
}
 
Matrix scalarMultiply(const double& value) {
if ( column_count != 1 ) {
throw std::invalid_argument("Incompatible matrix dimension.");
}
 
Matrix result(row_count, column_count);
for ( int32_t i = 0; i < row_count; ++i ) {
result.data[i][0] = data[i][0] * value;
}
return result;
}
 
Matrix unit() {
if ( column_count != 1 ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
const double the_magnitude = magnitude();
Matrix result(row_count, column_count);
for ( int32_t i = 0; i < row_count; ++i ) {
result.data[i][0] = data[i][0] / the_magnitude;
}
return result;
}
 
double magnitude() {
if ( column_count != 1 ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
double norm = 0.0;
for ( int32_t i = 0; i < row_count; ++i ) {
norm += data[i][0] * data[i][0];
}
return std::sqrt(norm);
}
 
int32_t size() {
if ( column_count != 1 ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
return row_count;
}
 
void display(const std::string& title) {
std::cout << title << std::endl;
for ( int32_t i = 0; i < row_count; ++i ) {
for ( int32_t j = 0; j < column_count; ++j ) {
std::cout << std::setw(9) << std::fixed << std::setprecision(4) << data[i][j];
}
std::cout << std::endl;
}
std::cout << std::endl;
}
 
double get_entry(const int32_t& row, const int32_t& col) {
return data[row][col];
}
 
void set_entry(const int32_t& row, const int32_t& col, const double& value) {
data[row][col] = value;
}
 
int32_t get_row_count() {
return row_count;
}
 
int32_t get_column_count() {
return column_count;
}
 
private:
void initialise() {
row_count = data.size();
column_count = data[0].size();
}
 
int32_t row_count;
int32_t column_count;
std::vector<std::vector<double>> data;
};
 
typedef std::pair<Matrix, Matrix> matrix_pair;
 
Matrix householder_factor(Matrix vector) {
if ( vector.get_column_count() != 1 ) {
throw std::invalid_argument("Incompatible matrix dimensions.");
}
 
const int32_t size = vector.size();
Matrix result(size, size);
for ( int32_t i = 0; i < size; ++i ) {
for ( int32_t j = 0; j < size; ++j ) {
result.set_entry(i, j, -2 * vector.get_entry(i, 0) * vector.get_entry(j, 0));
}
}
 
for ( int32_t i = 0; i < size; ++i ) {
result.set_entry(i, i, result.get_entry(i, i) + 1.0);
}
return result;
}
 
matrix_pair householder(Matrix matrix) {
const int32_t row_count = matrix.get_row_count();
const int32_t column_count = matrix.get_column_count();
std::vector<Matrix> versions_of_Q;
Matrix z(matrix);
Matrix z1(row_count, column_count);
 
for ( int32_t k = 0; k < column_count && k < row_count - 1; ++k ) {
Matrix vectorE(row_count, 1);
z1 = z.minor(k);
Matrix vectorX = z1.column(k);
double magnitudeX = vectorX.magnitude();
if ( matrix.get_entry(k, k) > 0 ) {
magnitudeX = -magnitudeX;
}
 
for ( int32_t i = 0; i < vectorE.size(); ++i ) {
vectorE.set_entry(i, 0, ( i == k ) ? 1 : 0);
}
vectorE = vectorE.scalarMultiply(magnitudeX).add(vectorX).unit();
versions_of_Q.emplace_back(householder_factor(vectorE));
z = versions_of_Q[k].multiply(z1);
}
 
Matrix Q = versions_of_Q[0];
for ( int32_t i = 1; i < column_count && i < row_count - 1; ++i ) {
Q = versions_of_Q[i].multiply(Q);
}
 
Matrix R = Q.multiply(matrix);
Q = Q.transpose();
return matrix_pair(R, Q);
}
 
Matrix solve_upper_triangular(Matrix r, Matrix b) {
const int32_t column_count = r.get_column_count();
Matrix result(column_count, 1);
 
for ( int32_t k = column_count - 1; k >= 0; --k ) {
double total = 0.0;
for ( int32_t j = k + 1; j < column_count; ++j ) {
total += r.get_entry(k, j) * result.get_entry(j, 0);
}
result.set_entry(k, 0, ( b.get_entry(k, 0) - total ) / r.get_entry(k, k));
}
return result;
}
 
Matrix least_squares(Matrix vandermonde, Matrix b) {
matrix_pair pair = householder(vandermonde);
return solve_upper_triangular(pair.first, pair.second.transpose().multiply(b));
}
 
Matrix fit_polynomial(Matrix x, Matrix y, const int32_t& polynomial_degree) {
Matrix vandermonde(x.get_column_count(), polynomial_degree + 1);
for ( int32_t i = 0; i < x.get_column_count(); ++i ) {
for ( int32_t j = 0; j < polynomial_degree + 1; ++j ) {
vandermonde.set_entry(i, j, std::pow(x.get_entry(0, i), j));
}
}
return least_squares(vandermonde, y.transpose());
}
 
int main() {
const std::vector<std::vector<double>> data = { { 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(data);
A.display("Initial matrix A:");
 
matrix_pair pair = householder(A);
Matrix Q = pair.second;
Matrix R = pair.first;
 
Q.display("Matrix Q:");
R.display("Matrix R:");
 
Matrix result = Q.multiply(R);
result.display("Matrix Q * R:");
 
// Task 2
Matrix x( std::vector<std::vector<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( std::vector<std::vector<double>>{
{ 1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0 } } );
 
result = fit_polynomial(x, y, 2);
result.display("Result of fitting polynomial:");
}
</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>
 
Line 2,950 ⟶ 3,269:
0,000 -175,000 70,000
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(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).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();
Matrix result = new Matrix(size, size);
for ( int i = 0; i < size; i++ ) {
for ( int j = 0; j < size; j++ ) {
result.setEntry(i, j, -2 * aVector.getEntry(i, 0) * aVector.getEntry(j, 0));
}
}
for ( int i = 0; i < size; i++ ) {
result.setEntry(i, i, result.getEntry(i, i) + 1.0);
}
return result;
}
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 aRowCount, int aColumnCount) {
this( new double[aRowCount][aColumnCount] );
}
 
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 aIndex) {
Matrix result = new Matrix(rowCount, 1);
for ( int i = 0; i < rowCount; i++ ) {
result.setEntry(i, 0, data[i][aIndex]);
}
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 < rowCount; 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 < rowCount; 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 rowCount;
}
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 aColumn) {
return data[aRow][aColumn];
}
public void setEntry(int aRow, int aColumn, double aValue) {
data[aRow][aColumn] = 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}}==
Line 4,804 ⟶ 5,442:
{{libheader|Wren-matrix}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="ecmascriptwren">import "./matrix" for Matrix
import "./fmt" for Fmt
 
var minor = Fn.new { |x, d|
PureFox
9,482

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/342343...350891"