September 2017 Update
This commit is contained in:
Ingy döt Net
parent
bba7bfd280
commit
ba8067c3b7
14570 changed files with 153136 additions and 63871 deletions
395
Task/QR-decomposition/C++/qr-decomposition.cpp
Normal file
395
Task/QR-decomposition/C++/qr-decomposition.cpp
Normal file
|
|
@ -0,0 +1,395 @@
|
|||
/*
|
||||
* g++ -O3 -Wall --std=c++11 qr_standalone.cpp -o qr_standalone
|
||||
*/
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <cstring> // for memset
|
||||
#include <limits>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
|
||||
#include <math.h>
|
||||
|
||||
class Vector;
|
||||
|
||||
class Matrix {
|
||||
|
||||
public:
|
||||
// default constructor (don't allocate)
|
||||
Matrix() : m(0), n(0), data(nullptr) {}
|
||||
|
||||
// constructor with memory allocation, initialized to zero
|
||||
Matrix(int m_, int n_) : Matrix() {
|
||||
m = m_;
|
||||
n = n_;
|
||||
allocate(m_,n_);
|
||||
}
|
||||
|
||||
// copy constructor
|
||||
Matrix(const Matrix& mat) : Matrix(mat.m,mat.n) {
|
||||
|
||||
for (int i = 0; i < m; i++)
|
||||
for (int j = 0; j < n; j++)
|
||||
(*this)(i,j) = mat(i,j);
|
||||
}
|
||||
|
||||
// constructor from array
|
||||
template<int rows, int cols>
|
||||
Matrix(double (&a)[rows][cols]) : Matrix(rows,cols) {
|
||||
|
||||
for (int i = 0; i < m; i++)
|
||||
for (int j = 0; j < n; j++)
|
||||
(*this)(i,j) = a[i][j];
|
||||
}
|
||||
|
||||
// destructor
|
||||
~Matrix() {
|
||||
deallocate();
|
||||
}
|
||||
|
||||
|
||||
// access data operators
|
||||
double& operator() (int i, int j) {
|
||||
return data[i+m*j]; }
|
||||
double operator() (int i, int j) const {
|
||||
return data[i+m*j]; }
|
||||
|
||||
// operator assignment
|
||||
Matrix& operator=(const Matrix& source) {
|
||||
|
||||
// self-assignment check
|
||||
if (this != &source) {
|
||||
if ( (m*n) != (source.m * source.n) ) { // storage cannot be reused
|
||||
allocate(source.m,source.n); // re-allocate storage
|
||||
}
|
||||
// storage can be used, copy data
|
||||
std::copy(source.data, source.data + source.m*source.n, data);
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// compute minor
|
||||
void compute_minor(const Matrix& mat, int d) {
|
||||
|
||||
allocate(mat.m, mat.n);
|
||||
|
||||
for (int i = 0; i < d; i++)
|
||||
(*this)(i,i) = 1.0;
|
||||
for (int i = d; i < mat.m; i++)
|
||||
for (int j = d; j < mat.n; j++)
|
||||
(*this)(i,j) = mat(i,j);
|
||||
|
||||
}
|
||||
|
||||
// Matrix multiplication
|
||||
// c = a * b
|
||||
// c will be re-allocated here
|
||||
void mult(const Matrix& a, const Matrix& b) {
|
||||
|
||||
if (a.n != b.m) {
|
||||
std::cerr << "Matrix multiplication not possible, sizes don't match !\n";
|
||||
return;
|
||||
}
|
||||
|
||||
// reallocate ourself if necessary i.e. current Matrix has not valid sizes
|
||||
if (a.m != m or b.n != n)
|
||||
allocate(a.m, b.n);
|
||||
|
||||
memset(data,0,m*n*sizeof(double));
|
||||
|
||||
for (int i = 0; i < a.m; i++)
|
||||
for (int j = 0; j < b.n; j++)
|
||||
for (int k = 0; k < a.n; k++)
|
||||
(*this)(i,j) += a(i,k) * b(k,j);
|
||||
|
||||
}
|
||||
|
||||
void transpose() {
|
||||
for (int i = 0; i < m; i++) {
|
||||
for (int j = 0; j < i; j++) {
|
||||
double t = (*this)(i,j);
|
||||
(*this)(i,j) = (*this)(j,i);
|
||||
(*this)(j,i) = t;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// take c-th column of m, put in v
|
||||
void extract_column(Vector& v, int c);
|
||||
|
||||
// memory allocation
|
||||
void allocate(int m_, int n_) {
|
||||
|
||||
// if already allocated, memory is freed
|
||||
deallocate();
|
||||
|
||||
// new sizes
|
||||
m = m_;
|
||||
n = n_;
|
||||
|
||||
data = new double[m_*n_];
|
||||
memset(data,0,m_*n_*sizeof(double));
|
||||
|
||||
} // allocate
|
||||
|
||||
// memory free
|
||||
void deallocate() {
|
||||
|
||||
if (data)
|
||||
delete[] data;
|
||||
|
||||
data = nullptr;
|
||||
|
||||
}
|
||||
|
||||
int m, n;
|
||||
|
||||
private:
|
||||
double* data;
|
||||
|
||||
}; // struct Matrix
|
||||
|
||||
// column vector
|
||||
class Vector {
|
||||
|
||||
public:
|
||||
// default constructor (don't allocate)
|
||||
Vector() : size(0), data(nullptr) {}
|
||||
|
||||
// constructor with memory allocation, initialized to zero
|
||||
Vector(int size_) : Vector() {
|
||||
size = size_;
|
||||
allocate(size_);
|
||||
}
|
||||
|
||||
// destructor
|
||||
~Vector() {
|
||||
deallocate();
|
||||
}
|
||||
|
||||
// access data operators
|
||||
double& operator() (int i) {
|
||||
return data[i]; }
|
||||
double operator() (int i) const {
|
||||
return data[i]; }
|
||||
|
||||
// operator assignment
|
||||
Vector& operator=(const Vector& source) {
|
||||
|
||||
// self-assignment check
|
||||
if (this != &source) {
|
||||
if ( size != (source.size) ) { // storage cannot be reused
|
||||
allocate(source.size); // re-allocate storage
|
||||
}
|
||||
// storage can be used, copy data
|
||||
std::copy(source.data, source.data + source.size, data);
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
// memory allocation
|
||||
void allocate(int size_) {
|
||||
|
||||
deallocate();
|
||||
|
||||
// new sizes
|
||||
size = size_;
|
||||
|
||||
data = new double[size_];
|
||||
memset(data,0,size_*sizeof(double));
|
||||
|
||||
} // allocate
|
||||
|
||||
// memory free
|
||||
void deallocate() {
|
||||
|
||||
if (data)
|
||||
delete[] data;
|
||||
|
||||
data = nullptr;
|
||||
|
||||
}
|
||||
|
||||
// ||x||
|
||||
double norm() {
|
||||
double sum = 0;
|
||||
for (int i = 0; i < size; i++) sum += (*this)(i) * (*this)(i);
|
||||
return sqrt(sum);
|
||||
}
|
||||
|
||||
// divide data by factor
|
||||
void rescale(double factor) {
|
||||
for (int i = 0; i < size; i++) (*this)(i) /= factor;
|
||||
}
|
||||
|
||||
void rescale_unit() {
|
||||
double factor = norm();
|
||||
rescale(factor);
|
||||
}
|
||||
|
||||
int size;
|
||||
|
||||
private:
|
||||
double* data;
|
||||
|
||||
}; // class Vector
|
||||
|
||||
// c = a + b * s
|
||||
void vmadd(const Vector& a, const Vector& b, double s, Vector& c)
|
||||
{
|
||||
if (c.size != a.size or c.size != b.size) {
|
||||
std::cerr << "[vmadd]: vector sizes don't match\n";
|
||||
return;
|
||||
}
|
||||
|
||||
for (int i = 0; i < c.size; i++)
|
||||
c(i) = a(i) + s * b(i);
|
||||
}
|
||||
|
||||
// mat = I - 2*v*v^T
|
||||
// !!! m is allocated here !!!
|
||||
void compute_householder_factor(Matrix& mat, const Vector& v)
|
||||
{
|
||||
|
||||
int n = v.size;
|
||||
mat.allocate(n,n);
|
||||
for (int i = 0; i < n; i++)
|
||||
for (int j = 0; j < n; j++)
|
||||
mat(i,j) = -2 * v(i) * v(j);
|
||||
for (int i = 0; i < n; i++)
|
||||
mat(i,i) += 1;
|
||||
}
|
||||
|
||||
// take c-th column of a matrix, put results in Vector v
|
||||
void Matrix::extract_column(Vector& v, int c) {
|
||||
if (m != v.size) {
|
||||
std::cerr << "[Matrix::extract_column]: Matrix and Vector sizes don't match\n";
|
||||
return;
|
||||
}
|
||||
|
||||
for (int i = 0; i < m; i++)
|
||||
v(i) = (*this)(i,c);
|
||||
}
|
||||
|
||||
void matrix_show(const Matrix& m, const std::string& str="")
|
||||
{
|
||||
std::cout << str << "\n";
|
||||
for(int i = 0; i < m.m; i++) {
|
||||
for (int j = 0; j < m.n; j++) {
|
||||
printf(" %8.3f", m(i,j));
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
|
||||
// L2-norm ||A-B||^2
|
||||
double matrix_compare(const Matrix& A, const Matrix& B) {
|
||||
// matrices must have same size
|
||||
if (A.m != B.m or A.n != B.n)
|
||||
return std::numeric_limits<double>::max();
|
||||
|
||||
double res=0;
|
||||
for(int i = 0; i < A.m; i++) {
|
||||
for (int j = 0; j < A.n; j++) {
|
||||
res += (A(i,j)-B(i,j)) * (A(i,j)-B(i,j));
|
||||
}
|
||||
}
|
||||
|
||||
res /= A.m*A.n;
|
||||
return res;
|
||||
}
|
||||
|
||||
void householder(Matrix& mat,
|
||||
Matrix& R,
|
||||
Matrix& Q)
|
||||
{
|
||||
|
||||
int m = mat.m;
|
||||
int n = mat.n;
|
||||
|
||||
// array of factor Q1, Q2, ... Qm
|
||||
std::vector<Matrix> qv(m);
|
||||
|
||||
// temp array
|
||||
Matrix z(mat);
|
||||
Matrix z1;
|
||||
|
||||
for (int k = 0; k < n && k < m - 1; k++) {
|
||||
|
||||
Vector e(m), x(m);
|
||||
double a;
|
||||
|
||||
// compute minor
|
||||
z1.compute_minor(z, k);
|
||||
|
||||
// extract k-th column into x
|
||||
z1.extract_column(x, k);
|
||||
|
||||
a = x.norm();
|
||||
if (mat(k,k) > 0) a = -a;
|
||||
|
||||
for (int i = 0; i < e.size; i++)
|
||||
e(i) = (i == k) ? 1 : 0;
|
||||
|
||||
// e = x + a*e
|
||||
vmadd(x, e, a, e);
|
||||
|
||||
// e = e / ||e||
|
||||
e.rescale_unit();
|
||||
|
||||
// qv[k] = I - 2 *e*e^T
|
||||
compute_householder_factor(qv[k], e);
|
||||
|
||||
// z = qv[k] * z1
|
||||
z.mult(qv[k], z1);
|
||||
|
||||
}
|
||||
|
||||
Q = qv[0];
|
||||
|
||||
// after this loop, we will obtain Q (up to a transpose operation)
|
||||
for (int i = 1; i < n && i < m - 1; i++) {
|
||||
|
||||
z1.mult(qv[i], Q);
|
||||
Q = z1;
|
||||
|
||||
}
|
||||
|
||||
R.mult(Q, mat);
|
||||
Q.transpose();
|
||||
}
|
||||
|
||||
double in[][3] = {
|
||||
{ 12, -51, 4},
|
||||
{ 6, 167, -68},
|
||||
{ -4, 24, -41},
|
||||
{ -1, 1, 0},
|
||||
{ 2, 0, 3},
|
||||
};
|
||||
|
||||
int main()
|
||||
{
|
||||
Matrix A(in);
|
||||
Matrix Q, R;
|
||||
|
||||
matrix_show(A,"A");
|
||||
|
||||
// compute QR decompostion
|
||||
householder(A, R, Q);
|
||||
|
||||
matrix_show(Q,"Q");
|
||||
matrix_show(R,"R");
|
||||
|
||||
// compare Q*R to the original matrix A
|
||||
Matrix A_check;
|
||||
A_check.mult(Q, R);
|
||||
|
||||
// compute L2 norm ||A-A_check||^2
|
||||
double l2 = matrix_compare(A,A_check);
|
||||
|
||||
// display Q*R
|
||||
matrix_show(A_check, l2 < 1e-12 ? "A == Q * R ? yes" : "A == Q * R ? no");
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -5,7 +5,7 @@
|
|||
typedef struct {
|
||||
int m, n;
|
||||
double ** v;
|
||||
} mat_t, mat;
|
||||
} mat_t, *mat;
|
||||
|
||||
mat matrix_new(int m, int n)
|
||||
{
|
||||
|
|
@ -37,7 +37,7 @@ void matrix_transpose(mat m)
|
|||
}
|
||||
}
|
||||
|
||||
mat matrix_copy(int n;double a[][n], int m, int n)
|
||||
mat matrix_copy(int n, double a[][n], int m)
|
||||
{
|
||||
mat x = matrix_new(m, n);
|
||||
for (int i = 0; i < m; i++)
|
||||
|
|
@ -174,7 +174,7 @@ double in[][3] = {
|
|||
int main()
|
||||
{
|
||||
mat R, Q;
|
||||
mat x = matrix_copy(in, 5, 3);
|
||||
mat x = matrix_copy(3, in, 5);
|
||||
householder(x, &R, &Q);
|
||||
|
||||
puts("Q"); matrix_show(Q);
|
||||
|
|
|
|||
89
Task/QR-decomposition/SequenceL/qr-decomposition.sequencel
Normal file
89
Task/QR-decomposition/SequenceL/qr-decomposition.sequencel
Normal file
|
|
@ -0,0 +1,89 @@
|
|||
import <Utilities/Math.sl>;
|
||||
import <Utilities/Sequence.sl>;
|
||||
import <Utilities/Conversion.sl>;
|
||||
|
||||
main :=
|
||||
let
|
||||
qrTest := [[12.0, -51.0, 4.0],
|
||||
[ 6.0, 167.0, -68.0],
|
||||
[-4.0, 24.0, -41.0]];
|
||||
|
||||
qrResult := qr(qrTest);
|
||||
|
||||
x := 1.0*(0 ... 10);
|
||||
y := 1.0*[1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
|
||||
|
||||
regResult := polyfit(x, y, 2);
|
||||
in
|
||||
"q:\n" ++ delimit(delimit(floatToString(qrResult[1], 6), ','), '\n') ++ "\n\n" ++
|
||||
"r:\n" ++ delimit(delimit(floatToString(qrResult[2], 1), ','), '\n') ++ "\n\n" ++
|
||||
"polyfit:\n" ++ "[" ++ delimit(floatToString(regResult, 1), ',') ++ "]";
|
||||
|
||||
//---Polynomial Regression---
|
||||
|
||||
polyfit(x(1), y(1), n) :=
|
||||
let
|
||||
a[j] := x ^ j foreach j within 0 ... n;
|
||||
in
|
||||
lsqr(transpose(a), transpose([y]));
|
||||
|
||||
lsqr(a(2), b(2)) :=
|
||||
let
|
||||
qrDecomp := qr(a);
|
||||
prod := mm(transpose(qrDecomp[1]), b);
|
||||
in
|
||||
solveUT(qrDecomp[2], prod);
|
||||
|
||||
solveUT(r(2), b(2)) :=
|
||||
let
|
||||
n := size(r[1]);
|
||||
in
|
||||
solveUTHelper(r, b, n, duplicate(0.0, n));
|
||||
|
||||
solveUTHelper(r(2), b(2), k, x(1)) :=
|
||||
let
|
||||
n := size(r[1]);
|
||||
newX := setElementAt(x, k, (b[k][1] - sum(r[k][(k+1) ... n] * x[(k+1) ... n])) / r[k][k]);
|
||||
in
|
||||
x when k <= 0
|
||||
else
|
||||
solveUTHelper(r, b, k - 1, newX);
|
||||
|
||||
//---QR Decomposition---
|
||||
|
||||
qr(A(2)) := qrHelper(A, id(size(A)), 1);
|
||||
|
||||
qrHelper(A(2), Q(2), i) :=
|
||||
let
|
||||
m := size(A);
|
||||
n := size(A[1]);
|
||||
|
||||
householder := makeHouseholder(A[i ... m, i]);
|
||||
|
||||
H[j,k] :=
|
||||
householder[j - i + 1][k - i + 1] when j >= i and k >= i
|
||||
else
|
||||
1.0 when j = k else 0.0
|
||||
foreach j within 1 ... m,
|
||||
k within 1 ... m;
|
||||
in
|
||||
[Q,A] when i > (n - 1 when m = n else n)
|
||||
else
|
||||
qrHelper(mm(H, A), mm(Q, H), i + 1);
|
||||
|
||||
|
||||
makeHouseholder(a(1)) :=
|
||||
let
|
||||
v := [1.0] ++ tail(a / (a[1] + sqrt(sum(a ^ 2)) * sign(a[1])));
|
||||
|
||||
H := id(size(a)) - (2.0 / mm([v], transpose([v])))[1,1] * mm(transpose([v]), [v]);
|
||||
in
|
||||
H;
|
||||
|
||||
//---Utilities---
|
||||
|
||||
id(n)[i,j] := 1.0 when i = j else 0.0
|
||||
foreach i within 1 ... n,
|
||||
j within 1 ... n;
|
||||
|
||||
mm(A(2), B(2))[i,j] := sum( A[i] * transpose(B)[j] );
|
||||
26
Task/QR-decomposition/Stata/qr-decomposition.stata
Normal file
26
Task/QR-decomposition/Stata/qr-decomposition.stata
Normal file
|
|
@ -0,0 +1,26 @@
|
|||
mata
|
||||
: qrd(a=(12,-51,4\6,167,-68\-4,24,-41),q=.,r=.)
|
||||
|
||||
: a
|
||||
1 2 3
|
||||
+-------------------+
|
||||
1 | 12 -51 4 |
|
||||
2 | 6 167 -68 |
|
||||
3 | -4 24 -41 |
|
||||
+-------------------+
|
||||
|
||||
: q
|
||||
1 2 3
|
||||
+----------------------------------------------+
|
||||
1 | -.8571428571 .3942857143 .3314285714 |
|
||||
2 | -.4285714286 -.9028571429 -.0342857143 |
|
||||
3 | .2857142857 -.1714285714 .9428571429 |
|
||||
+----------------------------------------------+
|
||||
|
||||
: r
|
||||
1 2 3
|
||||
+----------------------+
|
||||
1 | -14 -21 14 |
|
||||
2 | 0 -175 70 |
|
||||
3 | 0 0 -35 |
|
||||
+----------------------+
|
||||
8
Task/QR-decomposition/Zkl/qr-decomposition.zkl
Normal file
8
Task/QR-decomposition/Zkl/qr-decomposition.zkl
Normal file
|
|
@ -0,0 +1,8 @@
|
|||
var [const] GSL=Import("zklGSL"); // libGSL (GNU Scientific Library)
|
||||
A:=GSL.Matrix(3,3).set(12.0, -51.0, 4.0,
|
||||
6.0, 167.0, -68.0,
|
||||
4.0, 24.0, -41.0);
|
||||
Q,R:=A.QRDecomp();
|
||||
println("Q:\n",Q.format());
|
||||
println("R:\n",R.format());
|
||||
println("Q*R:\n",(Q*R).format());
|
||||
Loading…
Add table
Add a link
Reference in a new issue