Skip to content

Commit 5e7b86f

Browse files
ishii-norimiishii-norimi
authored
Add some iterative method of solve in Matrix class (#881)
1 parent 9dc949c commit 5e7b86f

File tree

2 files changed

+326
-0
lines changed

2 files changed

+326
-0
lines changed

‎lib/util/matrix.js‎

Lines changed: 149 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1814,6 +1814,23 @@ export default class Matrix {
18141814
return d
18151815
}
18161816

1817+
/**
1818+
* Returns a spectral radius.
1819+
* @returns {number} Spectral radius
1820+
*/
1821+
spectralRadius() {
1822+
if (!this.isSquare()) {
1823+
throw new MatrixException('Spectral radius only define square matrix.', this)
1824+
}
1825+
1826+
const eigVals = this.eigenValues()
1827+
let r = 0
1828+
for (let i = 0; i < eigVals.length; i++) {
1829+
r = Math.max(r, Math.abs(eigVals[i]))
1830+
}
1831+
return r
1832+
}
1833+
18171834
/**
18181835
* Multiply all elements by -1 in-place.
18191836
*/
@@ -3066,6 +3083,138 @@ export default class Matrix {
30663083
return x
30673084
}
30683085

3086+
/**
3087+
* Returns a solved value with Jacobi method.
3088+
* @param {Matrix} b Dependent variable values
3089+
* @param {number} [maxIteration] Maximum iteration
3090+
* @returns {Matrix} Solved matrix
3091+
*/
3092+
solveJacobi(b, maxIteration = 1.0e3) {
3093+
if (!this.isSquare()) {
3094+
throw new MatrixException('solveJacobi only define square matrix.', this)
3095+
}
3096+
if (this.rows !== b.rows) {
3097+
throw new MatrixException('b size is invalid.', [this, b])
3098+
}
3099+
const n = this.cols
3100+
const m = b.cols
3101+
3102+
let x = Matrix.zeros(n, m)
3103+
for (let t = 0; t < maxIteration; t++) {
3104+
const newx = Matrix.zeros(n, m)
3105+
for (let i = 0; i < n; i++) {
3106+
const ii = this.at(i, i)
3107+
for (let k = 0; k < m; k++) {
3108+
let v = b.at(i, k)
3109+
for (let j = 0; j < n; j++) {
3110+
if (i === j) continue
3111+
v -= x.at(j, k) * this.at(i, j)
3112+
}
3113+
newx.set(i, k, v / ii)
3114+
}
3115+
}
3116+
if (newx.some(v => isNaN(v))) {
3117+
throw new MatrixException('Can not calculate solved value.', this)
3118+
}
3119+
const e = Matrix.sub(x, newx).norm()
3120+
if (e < 1.0e-8) {
3121+
x = newx
3122+
break
3123+
}
3124+
x = newx
3125+
}
3126+
return x
3127+
}
3128+
3129+
/**
3130+
* Returns a solved value with Gauss-Seidel method.
3131+
* @param {Matrix} b Dependent variable values
3132+
* @param {number} [maxIteration] Maximum iteration
3133+
* @returns {Matrix} Solved matrix
3134+
*/
3135+
solveGaussSeidel(b, maxIteration = 1.0e3) {
3136+
if (!this.isSquare()) {
3137+
throw new MatrixException('solveGaussSeidel only define square matrix.', this)
3138+
}
3139+
if (this.rows !== b.rows) {
3140+
throw new MatrixException('b size is invalid.', [this, b])
3141+
}
3142+
const n = this.cols
3143+
const m = b.cols
3144+
3145+
let x = Matrix.zeros(n, m)
3146+
for (let t = 0; t < maxIteration; t++) {
3147+
let e = 0
3148+
for (let i = 0; i < n; i++) {
3149+
const ii = this.at(i, i)
3150+
for (let k = 0; k < m; k++) {
3151+
let v = b.at(i, k)
3152+
for (let j = 0; j < n; j++) {
3153+
if (i === j) continue
3154+
v -= x.at(j, k) * this.at(i, j)
3155+
}
3156+
v /= ii
3157+
e += (x.at(i, k) - v) ** 2
3158+
x.set(i, k, v)
3159+
}
3160+
}
3161+
if (x.some(v => isNaN(v))) {
3162+
throw new MatrixException('Can not calculate solved value.', this)
3163+
}
3164+
if (Math.sqrt(e) < 1.0e-8) {
3165+
break
3166+
}
3167+
}
3168+
return x
3169+
}
3170+
3171+
/**
3172+
* Returns a solved value with Successive Over-Relaxation method.
3173+
* @param {Matrix} b Dependent variable values
3174+
* @param {Matrix} w Relaxation factor
3175+
* @param {number} [maxIteration] Maximum iteration
3176+
* @returns {Matrix} Solved matrix
3177+
*/
3178+
solveSOR(b, w, maxIteration = 1.0e3) {
3179+
if (!this.isSquare()) {
3180+
throw new MatrixException('solveSOR only define square matrix.', this)
3181+
}
3182+
if (this.rows !== b.rows) {
3183+
throw new MatrixException('b size is invalid.', [this, b])
3184+
}
3185+
if (w <= 0 || 2 <= w) {
3186+
throw new MatrixException('w must be positive and less than 2.', [this, w])
3187+
}
3188+
const n = this.cols
3189+
const m = b.cols
3190+
3191+
let x = Matrix.zeros(n, m)
3192+
for (let t = 0; t < maxIteration; t++) {
3193+
let e = 0
3194+
for (let i = 0; i < n; i++) {
3195+
const ii = this.at(i, i)
3196+
for (let k = 0; k < m; k++) {
3197+
let v = b.at(i, k)
3198+
for (let j = 0; j < n; j++) {
3199+
if (i === j) continue
3200+
v -= x.at(j, k) * this.at(i, j)
3201+
}
3202+
v /= ii
3203+
const xik = x.at(i, k)
3204+
e += (xik - v) ** 2
3205+
x.set(i, k, (1 - w) * xik + w * v)
3206+
}
3207+
}
3208+
if (x.some(v => isNaN(v))) {
3209+
throw new MatrixException('Can not calculate solved value.', this)
3210+
}
3211+
if (Math.sqrt(e) < 1.0e-8) {
3212+
break
3213+
}
3214+
}
3215+
return x
3216+
}
3217+
30693218
/**
30703219
* Returns a bidiagonal matrix.
30713220
* @returns {Matrix} Bidiagonal matrix

‎tests/lib/util/matrix.test.js‎

Lines changed: 177 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3277,6 +3277,27 @@ describe('Matrix', () => {
32773277
})
32783278
})
32793279

3280+
describe('spectralRadius', () => {
3281+
test('default', () => {
3282+
const mat = Matrix.randn(5, 5).gram()
3283+
const r = mat.spectralRadius()
3284+
expect(r).toBeGreaterThan(0)
3285+
3286+
const pmat = mat.copy()
3287+
const mmat = mat.copy()
3288+
for (let k = 0; k < 5; k++) {
3289+
pmat.subAt(k, k, r)
3290+
mmat.addAt(k, k, r)
3291+
}
3292+
expect(pmat.det() * mmat.det()).toBeCloseTo(0)
3293+
})
3294+
3295+
test('fail', () => {
3296+
const mat = new Matrix(5, 4)
3297+
expect(() => mat.spectralRadius()).toThrow('Spectral radius only define square matrix.')
3298+
})
3299+
})
3300+
32803301
test.each([
32813302
['not', v => +!v],
32823303
['bitnot', v => ~v],
@@ -4946,6 +4967,162 @@ describe('Matrix', () => {
49464967
})
49474968
})
49484969

4970+
describe('solveJacobi', () => {
4971+
test.each([0, 1, 5])('success %i', n => {
4972+
const x = Matrix.randn(n, n)
4973+
x.add(Matrix.eye(n, n, 100))
4974+
const b = Matrix.randn(n, 1)
4975+
4976+
const a = x.solveJacobi(b)
4977+
4978+
const t = x.dot(a)
4979+
for (let i = 0; i < b.rows; i++) {
4980+
for (let j = 0; j < b.cols; j++) {
4981+
expect(t.at(i, j)).toBeCloseTo(b.at(i, j))
4982+
}
4983+
}
4984+
})
4985+
4986+
test.each([0, 1, 5])('excessive columns (%i)', n => {
4987+
const x = Matrix.randn(n, n)
4988+
x.add(Matrix.eye(n, n, 100))
4989+
const b = Matrix.randn(n, 1 + Math.floor(Math.random() * 10))
4990+
4991+
const a = x.solveJacobi(b)
4992+
4993+
const t = x.dot(a)
4994+
for (let i = 0; i < b.rows; i++) {
4995+
for (let j = 0; j < b.cols; j++) {
4996+
expect(t.at(i, j)).toBeCloseTo(b.at(i, j))
4997+
}
4998+
}
4999+
})
5000+
5001+
test('fail invalid columns', () => {
5002+
const a = Matrix.randn(10, 4)
5003+
const b = Matrix.randn(10, 1)
5004+
expect(() => a.solveJacobi(b)).toThrow('solveJacobi only define square matrix.')
5005+
})
5006+
5007+
test('fail invalid rows', () => {
5008+
const a = Matrix.randn(3, 3)
5009+
const b = Matrix.randn(4, 1)
5010+
expect(() => a.solveJacobi(b)).toThrow('b size is invalid.')
5011+
})
5012+
5013+
test('fail can not calculate', () => {
5014+
const a = Matrix.randn(4, 4)
5015+
const b = Matrix.randn(4, 1)
5016+
expect(() => a.solveJacobi(b)).toThrow('Can not calculate solved value.')
5017+
})
5018+
})
5019+
5020+
describe('solveGaussSeidel', () => {
5021+
test.each([0, 1, 5])('success %i', n => {
5022+
const x = Matrix.randn(n, n)
5023+
x.add(Matrix.eye(n, n, 100))
5024+
const b = Matrix.randn(n, 1)
5025+
5026+
const a = x.solveGaussSeidel(b)
5027+
5028+
const t = x.dot(a)
5029+
for (let i = 0; i < b.rows; i++) {
5030+
for (let j = 0; j < b.cols; j++) {
5031+
expect(t.at(i, j)).toBeCloseTo(b.at(i, j))
5032+
}
5033+
}
5034+
})
5035+
5036+
test.each([0, 1, 5])('excessive columns (%i)', n => {
5037+
const x = Matrix.randn(n, n)
5038+
x.add(Matrix.eye(n, n, 100))
5039+
const b = Matrix.randn(n, 1 + Math.floor(Math.random() * 10))
5040+
5041+
const a = x.solveGaussSeidel(b)
5042+
5043+
const t = x.dot(a)
5044+
for (let i = 0; i < b.rows; i++) {
5045+
for (let j = 0; j < b.cols; j++) {
5046+
expect(t.at(i, j)).toBeCloseTo(b.at(i, j))
5047+
}
5048+
}
5049+
})
5050+
5051+
test('fail invalid columns', () => {
5052+
const a = Matrix.randn(10, 4)
5053+
const b = Matrix.randn(10, 1)
5054+
expect(() => a.solveGaussSeidel(b)).toThrow('solveGaussSeidel only define square matrix.')
5055+
})
5056+
5057+
test('fail invalid rows', () => {
5058+
const a = Matrix.randn(3, 3)
5059+
const b = Matrix.randn(4, 1)
5060+
expect(() => a.solveGaussSeidel(b)).toThrow('b size is invalid.')
5061+
})
5062+
5063+
test('fail can not calculate', () => {
5064+
const a = Matrix.randn(4, 4)
5065+
const b = Matrix.randn(4, 1)
5066+
expect(() => a.solveGaussSeidel(b)).toThrow('Can not calculate solved value.')
5067+
})
5068+
})
5069+
5070+
describe('solveSOR', () => {
5071+
test.each([0, 1, 5])('success %i', n => {
5072+
const x = Matrix.randn(n, n)
5073+
x.add(Matrix.eye(n, n, 100))
5074+
const b = Matrix.randn(n, 1)
5075+
5076+
const a = x.solveSOR(b, 0.8)
5077+
5078+
const t = x.dot(a)
5079+
for (let i = 0; i < b.rows; i++) {
5080+
for (let j = 0; j < b.cols; j++) {
5081+
expect(t.at(i, j)).toBeCloseTo(b.at(i, j))
5082+
}
5083+
}
5084+
})
5085+
5086+
test.each([0, 1, 5])('excessive columns (%i)', n => {
5087+
const x = Matrix.randn(n, n)
5088+
x.add(Matrix.eye(n, n, 100))
5089+
const b = Matrix.randn(n, 1 + Math.floor(Math.random() * 10))
5090+
5091+
const a = x.solveSOR(b, 0.8)
5092+
5093+
const t = x.dot(a)
5094+
for (let i = 0; i < b.rows; i++) {
5095+
for (let j = 0; j < b.cols; j++) {
5096+
expect(t.at(i, j)).toBeCloseTo(b.at(i, j))
5097+
}
5098+
}
5099+
})
5100+
5101+
test('fail invalid columns', () => {
5102+
const a = Matrix.randn(10, 4)
5103+
const b = Matrix.randn(10, 1)
5104+
expect(() => a.solveSOR(b, 0.8)).toThrow('solveSOR only define square matrix.')
5105+
})
5106+
5107+
test('fail invalid rows', () => {
5108+
const a = Matrix.randn(3, 3)
5109+
const b = Matrix.randn(4, 1)
5110+
expect(() => a.solveSOR(b, 0.8)).toThrow('b size is invalid.')
5111+
})
5112+
5113+
test.each([-1, 0, 2, 3])('fail invalid w', w => {
5114+
const a = Matrix.randn(4, 4)
5115+
const b = Matrix.randn(4, 1)
5116+
expect(() => a.solveSOR(b, w)).toThrow('w must be positive and less than 2.')
5117+
})
5118+
5119+
test('fail can not calculate', () => {
5120+
const a = Matrix.randn(4, 4)
5121+
const b = Matrix.randn(4, 1)
5122+
expect(() => a.solveSOR(b, 0.8)).toThrow('Can not calculate solved value.')
5123+
})
5124+
})
5125+
49495126
describe('bidiag', () => {
49505127
test.each([
49515128
[2, 3],

0 commit comments

Comments
 (0)