import { neumair_sum } from "../numerical/index.js";
import { simultaneous_poweriteration } from "../linear_algebra/index.js";
import { Randomizer } from "../util/index.js";
/**
* @class
* @alias Matrix
* @requires module:numerical/neumair_sum
*/
export class Matrix {
/**
* creates a new Matrix. Entries are stored in a Float64Array.
* @memberof module:matrix
* @param {number} rows - The amount of rows of the matrix.
* @param {number} cols - The amount of columns of the matrix.
* @param {(function|string|number)} value=0 - Can be a function with row and col as parameters, a number, or "zeros", "identity" or "I", or "center".
* - **function**: for each entry the function gets called with the parameters for the actual row and column.
* - **string**: allowed are
* - "zero", creates a zero matrix.
* - "identity" or "I", creates an identity matrix.
* - "center", creates an center matrix.
* - **number**: create a matrix filled with the given value.
* @example
*
* let A = new Matrix(10, 10, () => Math.random()); //creates a 10 times 10 random matrix.
* let B = new Matrix(3, 3, "I"); // creates a 3 times 3 identity matrix.
* @returns {Matrix} returns a {@link rows} times {@link cols} Matrix filled with {@link value}.
*/
constructor(rows = null, cols = null, value = null) {
this._rows = rows;
this._cols = cols;
this._data = null;
if (rows && cols) {
if (!value) {
this._data = new Float64Array(rows * cols);
return this;
}
if (typeof value === "function") {
this._data = new Float64Array(rows * cols);
for (let row = 0; row < rows; ++row) {
for (let col = 0; col < cols; ++col) {
this._data[row * cols + col] = value(row, col);
}
}
return this;
}
if (typeof value === "string") {
if (value === "zeros") {
return new Matrix(rows, cols, 0);
}
if (value === "identity" || value === "I") {
this._data = new Float64Array(rows * cols);
for (let row = 0; row < rows; ++row) {
this._data[row * cols + row] = 1;
}
return this;
}
if (value === "center" && rows == cols) {
this._data = new Float64Array(rows * cols);
value = (i, j) => (i === j ? 1 : 0) - 1 / rows;
for (let row = 0; row < rows; ++row) {
for (let col = 0; col < cols; ++col) {
this._data[row * cols + col] = value(row, col);
}
}
return this;
}
}
if (typeof value === "number") {
this._data = new Float64Array(rows * cols);
for (let row = 0; row < rows; ++row) {
for (let col = 0; col < cols; ++col) {
this._data[row * cols + col] = value;
}
}
return this;
}
}
return this;
}
/**
* Creates a Matrix out of {@link A}.
* @param {(Matrix|Array|Float64Array|number)} A - The matrix, array, or number, which should converted to a Matrix.
* @param {"row"|"col"|"diag"} [type = "row"] - If {@link A} is a Array or Float64Array, then type defines if it is a row- or a column vector.
* @returns {Matrix}
*
* @example
* let A = Matrix.from([[1, 0], [0, 1]]); //creates a two by two identity matrix.
* let S = Matrix.from([1, 2, 3], "diag"); // creates a 3 by 3 matrix with 1, 2, 3 on its diagonal. [[1, 0, 0], [0, 2, 0], [0, 0, 3]]
*/
static from(A, type = "row") {
if (A instanceof Matrix) {
return A.clone();
} else if (Matrix.isArray(A)) {
let m = A.length;
if (m === 0) throw new Error("Array is empty");
// 1d
if (!Matrix.isArray(A[0])) {
if (type === "row") {
return new Matrix(1, m, (_, j) => A[j]);
} else if (type === "col") {
return new Matrix(m, 1, (i) => A[i]);
} else if (type === "diag") {
return new Matrix(m, m, (i, j) => (i == j ? A[i] : 0));
} else {
throw new Error("1d array has NaN entries");
}
// 2d
} else {
let n = A[0].length;
for (let row = 0; row < m; ++row) {
if (A[row].length !== n) {
throw new Error("various array lengths");
}
}
return new Matrix(m, n, (i, j) => A[i][j]);
}
} else if (typeof A === "number") {
return new Matrix(1, 1, A);
} else {
throw new Error("error");
}
}
/**
* Returns the {@link row}<sup>th</sup> row from the Matrix.
* @param {Number} row
* @returns {Float64Array}
*/
row(row) {
const data = this.values;
const cols = this._cols;
return data.subarray(row * cols, (row + 1) * cols);
}
/**
* Returns an generator yielding each row of the Matrix.
* @yields {Float64Array}
*/
*iterate_rows() {
const cols = this._cols;
const rows = this._rows;
const data = this.values;
for (let row = 0; row < rows; ++row) {
yield data.subarray(row * cols, (row + 1) * cols);
}
}
/**
* Makes a {@link Matrix} object an iterable object.
* @yields {Float64Array}
*/
*[Symbol.iterator]() {
for (const row of this.iterate_rows()) {
yield row;
}
}
/**
* Sets the entries of {@link row}<sup>th</sup> row from the Matrix to the entries from {@link values}.
* @param {Number} row
* @param {Array} values
* @returns {Matrix}
*/
set_row(row, values) {
const cols = this._cols;
if (Matrix.isArray(values) && values.length === cols) {
const offset = row * cols;
for (let col = 0; col < cols; ++col) {
this.values[offset + col] = values[col];
}
} else if (values instanceof Matrix && values.shape[1] === cols && values.shape[0] === 1) {
const offset = row * cols;
for (let col = 0; col < cols; ++col) {
this.values[offset + col] = values._data[col];
}
} else {
throw new Error("Values not valid! Needs to be either an Array, a Float64Array, or a fitting Matrix!")
}
return this;
}
/**
* Swaps the rows {@link row1} and {@link row2} of the Matrix.
* @param {Number} row1
* @param {Number} row2
* @returns {Matrix}
*/
swap_rows(row1, row2) {
const cols = this._cols;
const data = this.values;
for (let i = row1 * cols, j = row2 * cols, col = 0; col < cols; ++col, ++i, ++j) {
const t = data[i];
data[i] = data[j];
data[j] = t;
}
}
/**
* Returns the {@link col}<sup>th</sup> column from the Matrix.
* @param {Number} col
* @returns {Array}
*/
col(col) {
const result_col = new Float64Array(this._rows);
for (let row = 0; row < this._rows; ++row) {
result_col[row] = this.values[row * this._cols + col];
}
return result_col;
}
/**
* Returns the {@link col}<sup>th</sup> entry from the {@link row}<sup>th</sup> row of the Matrix.
* @param {int} row
* @param {int} col
* @returns {float64}
*/
entry(row, col) {
return this.values[row * this._cols + col];
}
/**
* Sets the {@link col}<sup>th</sup> entry from the {@link row}<sup>th</sup> row of the Matrix to the given {@link value}.
* @param {int} row
* @param {int} col
* @param {float64} value
* @returns {Matrix}
*/
set_entry(row, col, value) {
this.values[row * this._cols + col] = value;
return this;
}
/**
* Adds a given {@link value} to the {@link col}<sup>th</sup> entry from the {@link row}<sup>th</sup> row of the Matrix.
* @param {int} row
* @param {int} col
* @param {float64} value
* @returns {Matrix}
*/
add_entry(row, col, value) {
this.values[row * this._cols + col] += value;
return this;
}
/**
* Subtracts a given {@link value} from the {@link col}<sup>th</sup> entry from the {@link row}<sup>th</sup> row of the Matrix.
* @param {int} row
* @param {int} col
* @param {float64} value
* @returns {Matrix}
*/
sub_entry(row, col, value) {
this.values[row * this._cols + col] -= value;
return this;
}
/**
* Returns a new transposed Matrix.
* @returns {Matrix}
*/
transpose() {
let B = new Matrix(this._cols, this._rows, (row, col) => this.entry(col, row));
return B;
}
/**
* Returns a new transposed Matrix. Short-form of {@function transpose}.
* @returns {Matrix}
*/
get T() {
return this.transpose();
}
/**
* Returns the inverse of the Matrix.
* @returns {Matrix}
*/
inverse() {
const rows = this._rows;
const cols = this._cols;
const A = this.clone();
const B = new Matrix(rows, cols, 'I');
// foreach column
for (let col = 0; col < cols; ++col) {
// Search for maximum in this column (pivot)
let max_idx = col;
let max_val = Math.abs(A.entry(col, col));
for (let row = col + 1; row < rows; ++row) {
const val = Math.abs(A.entry(row, col));
if (max_val < val) {
max_idx = row;
max_val = val;
}
}
if (max_val === 0) {
throw new Error('Cannot compute inverse of Matrix, determinant is zero');
}
// Swap maximum row with current row
if (max_idx !== col) {
A.swap_rows(col, max_idx);
B.swap_rows(col, max_idx);
}
// eliminate non-zero values on the other rows at column c
const A_col = A.row(col);
const B_col = B.row(col);
for (let row = 0; row < rows; ++row) {
if (row !== col) {
// eliminate value at column c and row r
const A_row = A.row(row);
const B_row = B.row(row);
if (A_row[col] !== 0) {
const f = A_row[col] / A_col[col];
// sub (f * row c) from row r to eliminate the value at column c
for (let s = col; s < cols; ++s) {
A_row[s] -= (f * A_col[s]);
}
for (let s = 0; s < cols; ++s) {
B_row[s] -= (f * B_col[s]);
}
}
} else {
// normalize value at Acc to 1 (diagonal):
// divide each value of row r=c by the value at Acc
const f = A_col[col];
for (let s = col; s < cols; ++s) {
A_col[s] /= f;
}
for (let s = 0; s < cols; ++s) {
B_col[s] /= f;
}
}
}
}
return B;
}
/**
* Returns the dot product. If {@link B} is an Array or Float64Array then an Array gets returned. If {@link B} is a Matrix then a Matrix gets returned.
* @param {(Matrix|Array|Float64Array)} B the right side
* @returns {(Matrix|Array)}
*/
dot(B) {
if (B instanceof Matrix) {
let A = this;
const [rows_A, cols_A] = A.shape;
const [rows_B, cols_B] = B.shape;
if (cols_A !== rows_B) {
throw new Error(`A.dot(B): A is a ${A.shape.join(" ⨯ ")}-Matrix, B is a ${B.shape.join(" ⨯ ")}-Matrix:
A has ${cols_A} cols and B ${rows_B} rows.
Must be equal!`);
}
const C = new Matrix(rows_A, cols_B, (row, col) => {
const A_i = A.row(row);
const B_val = B.values;
let sum = 0;
for (let i = 0, j = col; i < cols_A; ++i, j += cols_B) {
sum += A_i[i] * B_val[j];
}
return sum;
});
return C;
} else if (Matrix.isArray(B)) {
let rows = this._rows;
if (B.length !== rows) {
throw new Error(`A.dot(B): A has ${rows} cols and B has ${B.length} rows. Must be equal!`);
}
let C = new Array(rows);
for (let row = 0; row < rows; ++row) {
C[row] = neumair_sum(this.row(row).map((e) => e * B[row]));
}
return C;
} else {
throw new Error(`B must be Matrix or Array`);
}
}
/**
* Transposes the current matrix and returns the dot product with {@link B}.
* If {@link B} is an Array or Float64Array then an Array gets returned.
* If {@link B} is a Matrix then a Matrix gets returned.
* @param {(Matrix|Array|Float64Array)} B the right side
* @returns {(Matrix|Array)}
*/
transDot(B) {
if (B instanceof Matrix) {
let A = this;
const [cols_A, rows_A] = A.shape; // transpose matrix
const [rows_B, cols_B] = B.shape;
if (cols_A !== rows_B) {
throw new Error(`A.dot(B): A is a ${[rows_A, cols_A].join(" ⨯ ")}-Matrix, B is a ${B.shape.join(" ⨯ ")}-Matrix:
A has ${cols_A} cols and B ${rows_B} rows, which must be equal!`);
}
// let B = new Matrix(this._cols, this._rows, (row, col) => this.entry(col, row));
// this.values[row * this._cols + col];
const C = new Matrix(rows_A, cols_B, (row, col) => {
const A_val = A.values;
const B_val = B.values;
let sum = 0;
for (let i = 0, j = row, k = col; i < cols_A; ++i, j += rows_A, k += cols_B) {
sum += A_val[j] * B_val[k];
}
return sum;
});
return C;
} else if (Matrix.isArray(B)) {
let rows = this._cols;
if (B.length !== rows) {
throw new Error(`A.dot(B): A has ${rows} cols and B has ${B.length} rows. Must be equal!`);
}
let C = new Array(rows);
for (let row = 0; row < rows; ++row) {
C[row] = neumair_sum(this.col(row).map((e) => e * B[row]));
}
return C;
} else {
throw new Error(`B must be Matrix or Array`);
}
}
/**
* Returns the dot product with the transposed version of {@link B}.
* If {@link B} is an Array or Float64Array then an Array gets returned.
* If {@link B} is a Matrix then a Matrix gets returned.
* @param {(Matrix|Array|Float64Array)} B the right side
* @returns {(Matrix|Array)}
*/
dotTrans(B) {
if (B instanceof Matrix) {
let A = this;
const [rows_A, cols_A] = A.shape;
const [cols_B, rows_B] = B.shape;
if (cols_A !== rows_B) {
throw new Error(`A.dot(B): A is a ${A.shape.join(" ⨯ ")}-Matrix, B is a ${[rows_B, cols_B].join(" ⨯ ")}-Matrix:
A has ${cols_A} cols and B ${rows_B} rows, which must be equal!`);
}
const C = new Matrix(rows_A, cols_B, (row, col) => {
const A_i = A.row(row);
const B_i = B.row(col);
let sum = 0;
for (let i = 0; i < cols_A; ++i) {
sum += A_i[i] * B_i[i];
}
return sum;
});
return C;
} else if (Matrix.isArray(B)) {
let rows = this._rows;
if (B.length !== rows) {
throw new Error(`A.dot(B): A has ${rows} cols and B has ${B.length} rows. Must be equal!`);
}
let C = new Array(rows);
for (let row = 0; row < rows; ++row) {
C[row] = neumair_sum(this.row(row).map((e) => e * B[row]));
}
return C;
} else {
throw new Error(`B must be Matrix or Array`);
}
}
/**
* Computes the outer product from {@link this} and {@link B}.
* @param {Matrix} B
* @returns {Matrix}
*/
outer(B) {
let A = this;
let l = A._data.length;
let r = B._data.length;
if (l != r) return undefined;
let C = new Matrix();
C.shape = [
l,
l,
(i, j) => {
if (i <= j) {
return A._data[i] * B._data[j];
} else {
return C.entry(j, i);
}
},
];
return C;
}
/**
* Appends matrix {@link B} to the matrix.
* @param {Matrix} B - matrix to append.
* @param {"horizontal"|"vertical"|"diag"} [type = "horizontal"] - type of concatenation.
* @returns {Matrix}
* @example
*
* let A = Matrix.from([[1, 1], [1, 1]]); // 2 by 2 matrix filled with ones.
* let B = Matrix.from([[2, 2], [2, 2]]); // 2 by 2 matrix filled with twos.
*
* A.concat(B, "horizontal"); // 2 by 4 matrix. [[1, 1, 2, 2], [1, 1, 2, 2]]
* A.concat(B, "vertical"); // 4 by 2 matrix. [[1, 1], [1, 1], [2, 2], [2, 2]]
* A.concat(B, "diag"); // 4 by 4 matrix. [[1, 1, 0, 0], [1, 1, 0, 0], [0, 0, 2, 2], [0, 0, 2, 2]]
*/
concat(B, type = "horizontal") {
const A = this;
const [rows_A, cols_A] = A.shape;
const [rows_B, cols_B] = B.shape;
if (type == "horizontal") {
if (rows_A != rows_B) {
throw new Error(`A.concat(B, "horizontal"): A and B need same number of rows, A has ${rows_A} rows, B has ${rows_B} rows.`);
}
const X = new Matrix(rows_A, cols_A + cols_B, "zeros");
X.set_block(0, 0, A);
X.set_block(0, cols_A, B);
return X;
} else if (type == "vertical") {
if (cols_A != cols_B) {
throw new Error(`A.concat(B, "vertical"): A and B need same number of columns, A has ${cols_A} columns, B has ${cols_B} columns.`);
}
const X = new Matrix(rows_A + rows_B, cols_A, "zeros");
X.set_block(0, 0, A);
X.set_block(rows_A, 0, B);
return X;
} else if (type == "diag") {
const X = new Matrix(rows_A + rows_B, cols_A + cols_B, "zeros");
X.set_block(0, 0, A);
X.set_block(rows_A, cols_A, B);
return X;
} else {
throw new Error(`type must be "horizontal" or "vertical", but type is ${type}!`);
}
}
/**
* Writes the entries of B in A at an offset position given by {@link offset_row} and {@link offset_col}.
* @param {int} offset_row
* @param {int} offset_col
* @param {Matrix} B
* @returns {Matrix}
*/
set_block(offset_row, offset_col, B) {
const rows = Math.min(this._rows - offset_row, B.shape[0]);
const cols = Math.min(this._cols - offset_col, B.shape[1]);
for (let row = 0; row < rows; ++row) {
for (let col = 0; col < cols; ++col) {
this.set_entry(row + offset_row, col + offset_col, B.entry(row, col));
}
}
return this;
}
/**
* Extracts the entries from the {@link start_row}<sup>th</sup> row to the {@link end_row}<sup>th</sup> row, the {@link start_col}<sup>th</sup> column to the {@link end_col}<sup>th</sup> column of the matrix.
* If {@link end_row} or {@link end_col} is empty, the respective value is set to {@link this.rows} or {@link this.cols}.
* @param {Number} start_row
* @param {Number} start_col
* @param {Number} [end_row = null]
* @param {Number} [end_col = null]
* @returns {Matrix} Returns a end_row - start_row times end_col - start_col matrix, with respective entries from the matrix.
* @example
*
* let A = Matrix.from([[1, 2, 3], [4, 5, 6], [7, 8, 9]]); // a 3 by 3 matrix.
*
* A.get_block(1, 1); // [[5, 6], [8, 9]]
* A.get_block(0, 0, 1, 1); // [[1]]
* A.get_block(1, 1, 2, 2); // [[5]]
* A.get_block(0, 0, 2, 2); // [[1, 2], [4, 5]]
*/
get_block(start_row, start_col, end_row = null, end_col = null) {
const [rows, cols] = this.shape;
end_row = end_row ?? rows;
end_col = end_col ?? cols;
if (end_row <= start_row || end_col <= start_col) {
throw new Error(`
end_row must be greater than start_row, and
end_col must be greater than start_col, but
end_row = ${end_row}, start_row = ${start_row}, end_col = ${end_col}, and start_col = ${start_col}!`);
}
const X = new Matrix(end_row - start_row, end_col - start_col, "zeros");
for (let row = start_row, new_row = 0; row < end_row; ++row, ++new_row) {
for (let col = start_col, new_col = 0; col < end_col; ++col, ++new_col) {
X.set_entry(new_row, new_col, this.entry(row, col));
}
}
return X;
//return new Matrix(end_row - start_row, end_col - start_col, (i, j) => this.entry(i + start_row, j + start_col));
}
/**
* Returns a new array gathering entries defined by the indices given by argument.
* @param {Array<Number>} row_indices - Array consists of indices of rows for gathering entries of this matrix
* @param {Array<Number>} col_indices - Array consists of indices of cols for gathering entries of this matrix
* @returns {Matrix}
*/
gather(row_indices, col_indices) {
const N = row_indices.length;
const D = col_indices.length;
const R = new Matrix(N, D);
for (let i = 0; i < N; ++i) {
const row_index = row_indices[i];
for (let j = 0; j < N; ++j) {
const col_index = col_indices[j];
R.set_entry(i, j, this.entry(row_index, col_index));
}
}
return R;
}
/**
* Applies a function to each entry of the matrix.
* @private
* @param {Function} f function takes 2 parameters, the value of the actual entry and a value given by the function {@link v}. The result of {@link f} gets writen to the Matrix.
* @param {Function} v function takes 2 parameters for row and col, and returns a value witch should be applied to the colth entry of the rowth row of the matrix.
*/
_apply_array(f, v) {
const data = this.values;
const [rows, cols] = this.shape;
for (let i = 0, row = 0; row < rows; ++row) {
for (let col = 0; col < cols; ++col, ++i) {
data[i] = f(data[i], v(row, col));
}
}
return this;
}
_apply_rowwise_array(values, f) {
return this._apply_array(f, (_, j) => values[j]);
}
_apply_colwise_array(values, f) {
const data = this.values;
const [rows, cols] = this.shape;
for (let i = 0, row = 0; row < rows; ++row) {
const val = values[row];
for (let col = 0; col < cols; ++col, ++i) {
data[i] = f(data[i], val);
}
}
return this;
}
_apply(value, f) {
const data = this.values;
const [rows, cols] = this.shape;
if (value instanceof Matrix) {
const values = value.values;
const [value_rows, value_cols] = value.shape;
if (value_rows === 1) {
if (cols !== value_cols) {
throw new Error(`cols !== value_cols`);
}
for (let i = 0, row = 0; row < rows; ++row) {
for (let col = 0; col < cols; ++col, ++i) {
data[i] = f(data[i], values[col]);
}
}
} else if (value_cols === 1) {
if (rows !== value_rows) {
throw new Error(`rows !== value_rows`);
}
for (let i = 0, row = 0; row < rows; ++row) {
const v = values[row];
for (let col = 0; col < cols; ++col, ++i) {
data[i] = f(data[i], v);
}
}
} else if (rows == value_rows && cols == value_cols) {
for (let i = 0, n = rows * cols; i < n; ++i) {
data[i] = f(data[i], values[i]);
}
} else {
throw new Error(`error`);
}
} else if (Matrix.isArray(value)) {
if (value.length === rows) {
for (let i = 0, row = 0; row < rows; ++row) {
const v = value[row];
for (let col = 0; col < cols; ++col, ++i) {
data[i] = f(data[i], v);
}
}
} else if (value.length === cols) {
for (let i = 0, row = 0; row < rows; ++row) {
for (let col = 0; col < cols; ++col, ++i) {
data[i] = f(data[i], value[col]);
}
}
} else {
throw new Error(`error`);
}
} else { // scalar value
for (let i = 0, n = rows * cols; i < n; ++i) {
data[i] = f(data[i], value);
}
}
return this;
}
/**
* Clones the Matrix.
* @returns {Matrix}
*/
clone() {
let B = new Matrix();
B._rows = this._rows;
B._cols = this._cols;
B._data = this.values.slice(0);
return B;
}
/**
* Entrywise multiplication with {@link value}.
* @param {Matrix|Array|Number} value
* @param {Object} [options]
* @param {Boolean} [options.inline = false] - If true, applies multiplication to the element, otherwise it creates first a copy and applies the multiplication on the copy.
* @returns {Matrix}
* @example
*
* let A = Matrix.from([[1, 2], [3, 4]]); // a 2 by 2 matrix.
* let B = A.clone(); // B == A;
*
* A.mult(2); // [[2, 4], [6, 8]];
* A.mult(B); // [[1, 4], [9, 16]];
*/
mult(value, { inline = false } = {}) {
const A = inline ? this : this.clone();
return A._apply(value, (a, b) => a * b);
}
/**
* Entrywise division with {@link value}.
* @param {Matrix|Array|Number} value
* @param {Object} [options]
* @param {Boolean} [options.inline = false] - If true, applies division to the element, otherwise it creates first a copy and applies the division on the copy.
* @returns {Matrix}
* @example
*
* let A = Matrix.from([[1, 2], [3, 4]]); // a 2 by 2 matrix.
* let B = A.clone(); // B == A;
*
* A.divide(2); // [[0.5, 1], [1.5, 2]];
* A.divide(B); // [[1, 1], [1, 1]];
*/
divide(value, { inline = false } = {}) {
const A = inline ? this : this.clone();
return A._apply(value, (a, b) => a / b);
}
/**
* Entrywise addition with {@link value}.
* @param {Matrix|Array|Number} value
* @param {Object} [options]
* @param {Boolean} [options.inline = false] - If true, applies addition to the element, otherwise it creates first a copy and applies the addition on the copy.
* @returns {Matrix}
* @example
*
* let A = Matrix.from([[1, 2], [3, 4]]); // a 2 by 2 matrix.
* let B = A.clone(); // B == A;
*
* A.add(2); // [[3, 4], [5, 6]];
* A.add(B); // [[2, 4], [6, 8]];
*/
add(value, {inline = false} = {}) {
const A = inline ? this : this.clone();
return A._apply(value, (a, b) => a + b);
}
/**
* Entrywise subtraction with {@link value}.
* @param {Matrix|Array|Number} value
* @param {Object} [options]
* @param {Boolean} [options.inline = false] - If true, applies subtraction to the element, otherwise it creates first a copy and applies the subtraction on the copy.
* @returns {Matrix}
* @example
*
* let A = Matrix.from([[1, 2], [3, 4]]); // a 2 by 2 matrix.
* let B = A.clone(); // B == A;
*
* A.sub(2); // [[-1, 0], [1, 2]];
* A.sub(B); // [[0, 0], [0, 0]];
*/
sub(value, { inline = false } = {}) {
const A = inline ? this : this.clone();
return A._apply(value, (a, b) => a - b);
}
/**
* Returns the number of rows and columns of the Matrix.
* @returns {Array} An Array in the form [rows, columns].
*/
get shape() {
return [this._rows, this._cols];
}
/**
* Returns the matrix in the given shape with the given function which returns values for the entries of the matrix.
* @param {Array} parameter - takes an Array in the form [rows, cols, value], where rows and cols are the number of rows and columns of the matrix, and value is a function which takes two parameters (row and col) which has to return a value for the colth entry of the rowth row.
* @returns {Matrix}
*/
set shape([rows, cols, value = () => 0]) {
this._rows = rows;
this._cols = cols;
this._data = new Float64Array(rows * cols);
for (let i = 0, row = 0; row < rows; ++row) {
for (let col = 0; col < cols; ++col, ++i) {
this._data[i] = value(row, col);
}
}
return this;
}
/**
* Returns the Matrix as a Array of Float64Arrays.
* @returns {Array<Float64Array>}
*/
get to2dArray() {
const result = [];
for (const row of this.iterate_rows()) {
result.push(row);
}
return result;
}
/**
* Returns the Matrix as a Array of Arrays.
* @returns {Array<Array>}
*/
get asArray() {
const result = [];
for (const row of this.iterate_rows()) {
result.push(Array.from(row));
}
return result;
}
/**
* Returns the diagonal of the Matrix.
* @returns {Float64Array}
*/
get diag() {
const rows = this._rows;
const cols = this._cols;
const min_row_col = Math.min(rows, cols);
let result = new Float64Array(min_row_col);
for (let i = 0; i < min_row_col; ++i) {
result[i] = this.entry(i, i);
}
return result;
}
/**
* Returns the mean of all entries of the Matrix.
* @returns {Number}
*/
get mean() {
const sum = this.sum;
const n = this._rows * this._cols;
return sum / n;
}
/**
* Returns the sum oof all entries of the Matrix.
* @returns {Number}
*/
get sum() {
const data = this.values;
return neumair_sum(data);
}
/**
* Returns the entries of the Matrix.
* @returns {Float64Array}
*/
get values() {
const data = this._data;
return data;
}
/**
* Returns the mean of each row of the matrix.
* @returns {Float64Array}
*/
get meanRows() {
const data = this.values;
const rows = this._rows;
const cols = this._cols;
const result = Float64Array.from({ length: rows });
for (let i = 0, row = 0; row < rows; ++row) {
let sum = 0;
for (let col = 0; col < cols; ++col, ++i) {
sum += data[i];
}
result[row] = sum / cols;
}
return result;
}
/** Returns the mean of each column of the matrix.
* @returns {Float64Array}
*/
get meanCols() {
const data = this.values;
const rows = this._rows;
const cols = this._cols;
const result = Float64Array.from({ length: cols });
for (let col = 0; col < cols; ++col) {
let sum = 0;
for (let i = col, row = 0; row < rows; ++row, i += cols) {
sum += data[i];
}
result[col] = sum / rows;
}
return result;
}
/**
* Solves the equation {@link A}x = {@link b} using the conjugate gradient method. Returns the result x.
* @param {Matrix} A - Matrix
* @param {Matrix} b - Matrix
* @param {Randomizer} [randomizer=null]
* @param {Number} [tol=1e-3]
* @returns {Matrix}
*/
static solve_CG(A, b, randomizer, tol = 1e-3) {
if (randomizer === null) {
randomizer = new Randomizer();
}
const rows = A.shape[0];
const cols = b.shape[1];
let result = new Matrix(rows, 0);
for (let i = 0; i < cols; ++i) {
const b_i = Matrix.from(b.col(i)).T;
let x = new Matrix(rows, 1, () => randomizer.random);
let r = b_i.sub(A.dot(x));
let d = r.clone();
do {
const z = A.dot(d);
const alpha = r.transDot(r).entry(0, 0) / d.transDot(z).entry(0, 0);
x = x.add(d.mult(alpha));
const r_next = r.sub(z.mult(alpha));
const beta = r_next.transDot(r_next).entry(0, 0) / r.transDot(r).entry(0, 0);
d = r_next.add(d.mult(beta));
r = r_next;
} while (Math.abs(r.mean) > tol);
result = result.concat(x, "horizontal");
}
return result;
}
/**
* Solves the equation {@link A}x = {@link b}. Returns the result x.
* @param {Matrix} A - Matrix or LU Decomposition
* @param {Matrix} b - Matrix
* @returns {Matrix}
*/
static solve(A, b) {
let { L: L, U: U } = "L" in A && "U" in A ? A : Matrix.LU(A);
let rows = L.shape[0];
let x = b.clone();
// forward
for (let row = 0; row < rows; ++row) {
for (let col = 0; col < row - 1; ++col) {
x.sub_entry(0, row, L.entry(row, col) * x.entry(1, col));
}
x.set_entry(0, row, x.entry(0, row) / L.entry(row, row));
}
// backward
for (let row = rows - 1; row >= 0; --row) {
for (let col = rows - 1; col > row; --col) {
x.sub_entry(0, row, U.entry(row, col) * x.entry(0, col));
}
x.set_entry(0, row, x.entry(0, row) / U.entry(row, row));
}
return x;
}
/**
* {@link L}{@link U} decomposition of the Matrix {@link A}. Creates two matrices, so that the dot product LU equals A.
* @param {Matrix} A
* @returns {{L: Matrix, U: Matrix}} result - Returns the left triangle matrix {@link L} and the upper triangle matrix {@link U}.
*/
static LU(A) {
const rows = A.shape[0];
const L = new Matrix(rows, rows, "zeros");
const U = new Matrix(rows, rows, "identity");
for (let j = 0; j < rows; ++j) {
for (let i = j; i < rows; ++i) {
let sum = 0;
for (let k = 0; k < j; ++k) {
sum += L.entry(i, k) * U.entry(k, j);
}
L.set_entry(i, j, A.entry(i, j) - sum);
}
for (let i = j; i < rows; ++i) {
if (L.entry(j, j) === 0) {
return undefined;
}
let sum = 0;
for (let k = 0; k < j; ++k) {
sum += L.entry(j, k) * U.entry(k, i);
}
U.set_entry(j, i, (A.entry(j, i) - sum) / L.entry(j, j));
}
}
return { L: L, U: U };
}
/**
* Computes the determinante of {@link A}, by using the LU decomposition of {@link A}.
* @param {Matrix} A
* @returns {Number} det - Returns the determinate of the Matrix {@link A}.
*/
static det(A) {
const rows = A.shape[0];
const { L, U } = Matrix.LU(A);
const L_diag = L.diag;
const U_diag = U.diag;
let det = L_diag[0] * U_diag[0];
for (let row = 1; row < rows; ++row) {
det *= L_diag[row] * U_diag[row];
}
return det;
}
/**
* Computes the {@link k} components of the SVD decomposition of the matrix {@link M}
* @param {Matrix} M
* @param {int} [k=2]
* @returns {{U: Matrix, Sigma: Matrix, V: Matrix}}
*/
static SVD(M, k = 2) {
let MtM = M.transDot(M);
let MMt = M.dotTrans(M);
let { eigenvectors: V, eigenvalues: Sigma } = simultaneous_poweriteration(MtM, k);
let { eigenvectors: U } = simultaneous_poweriteration(MMt, k);
return { U: U, Sigma: Sigma.map((sigma) => Math.sqrt(sigma)), V: V };
//Algorithm 1a: Householder reduction to bidiagonal form:
/* const [m, n] = A.shape;
let U = new Matrix(m, n, (i, j) => i == j ? 1 : 0);
console.log(U.to2dArray)
let V = new Matrix(n, m, (i, j) => i == j ? 1 : 0);
console.log(V.to2dArray)
let B = Matrix.bidiagonal(A.clone(), U, V);
console.log(U,V,B)
return { U: U, "Sigma": B, V: V }; */
}
static isArray(A) {
return Array.isArray(A) || A instanceof Float64Array || A instanceof Float32Array;
}
}