jacobi_alg

jacobi.js

import Matrix2 from '../core/mat2.js'
import Matrix3, { transpose, multiply, } from '../core/mat3.js'

/**
 * 雅可比迭代法求解实对称矩阵的特征向量、特征值
 * @param {Matrix2 | Matrix3} matrix 二维数组，当前支持二维、三维迭代求解
 * @returns {Object} 返回一个对象 obj, obj.eigenVectors 是特征向量(二维数组)，obj.eigenValues 是特征值(一维数组)
 */
function jacobi(matrix) {
  const length = matrix.length === undefined ? matrix.dimensions : matrix.length
  if (length === 2)
    return jacobi2d(matrix)
  else if (length === 3)
    return jacobi3d(matrix)
  else
    throw Error('当前还不支持其他维度的实对称矩阵求解特征向量、特征值')
}

function jacobi2d(matrix) {
  const m = Matrix2.fromElementsArray(matrix.flat())
  let rotationAngle = 0
  /** 计算 Givens 旋转矩阵的旋转角 */
  if (m.m11 - m.m22 === 0) {
    rotationAngle = Math.PI / 4
  } else {
    const tan2xita = 2 * m.m12 / (m.m11 - m.m22)
    rotationAngle = Math.atan2(tan2xita, 1) / 2
  }

  /** 构造二阶 Givens 旋转矩阵 */
  const givensMatrix = Matrix2.fromElementsArray([
    Math.cos(rotationAngle), Math.sin(rotationAngle),
    -1 * Math.sin(rotationAngle), Math.cos(rotationAngle)
  ])
  /** 获取 Givens 旋转矩阵的逆矩阵 */
  const inverseGivensMatrix = Matrix2.inverse(givensMatrix)
  /** 求解得到旋转后的矩阵，即斜对角线特征值矩阵 */
  const eigenValuesMatrix = Matrix2.multiply(Matrix2.multiply(givensMatrix, m), inverseGivensMatrix)

  return {
    eigenValues: [eigenValuesMatrix.m11, eigenValuesMatrix.m22],
    eigenVectors: [
      [inverseGivensMatrix.m11, inverseGivensMatrix.m21],
      [inverseGivensMatrix.m12, inverseGivensMatrix.m22]
    ]
  }
}

function getGivensIndexes(m) {
  const m12 = m.m12, m13 = m.m13, m23 = m.m23
  let max = m12 > m13 ? m12 : m13
  max = max > m23 ? max : m23
  if (max === m12) {
    return [1, 2]
  } else if (max === m13) {
    return [1, 3]
  } else {
    return [2, 3]
  }
}

function createGivens(m, p, q, result) {
  if (result === undefined) {
    result = new Matrix3()
  }

  let rad = 0
  const m_pp = m[`m${p}${p}`]
  const m_qq = m[`m${q}${q}`]
  const m_pq = m[`m${p}${q}`]
  if (m_pp === m_qq) {
    rad = Math.PI / 4
  } else {
    rad = Math.atan2(2 * m_pq / (m_pp - m_qq), 1) / 2
  }

  // 例如
  // [
  //   mpp, mpq, 0
  //   mqp, mqq, 0
  //   0,   0,   1
  // ]

  result[`m${p}${p}`] = Math.cos(rad)
  result[`m${q}${q}`] = Math.cos(rad)
  result[`m${p}${q}`] = Math.sin(rad)
  result[`m${q}${p}`] = -Math.sin(rad)
  return result
}

function jacobiLoopCondition(m) {
  const topTrianglesTimes2 = (m.m12 ** 2 + m.m13 ** 2 + m.m23 ** 2) * 2
  return topTrianglesTimes2 < 0.1 ** 10
}

let givensLeft = new Matrix3()
let givensRight = new Matrix3()
let times = 0
function jacobi3d(m, maxTime = 10) {
  times += 1
  /* step1. 寻找行列号 p q */
  const [p, q] = getGivensIndexes(m)
  /* step2. 构造 Givens 矩阵，及其逆矩阵（即转置） */
  const tempLeft = createGivens(m, p, q, new Matrix3())
  const tempRight = transpose(tempLeft, new Matrix3())

  /* step3. 累乘旧矩阵 */
  multiply(tempLeft, givensLeft, givensLeft)
  multiply(givensRight, tempRight, givensRight)

  /* step4. Givens 旋转 */
  const leftMultiply = multiply(tempLeft, m, new Matrix3())
  const tempMatrix = multiply(leftMultiply, tempRight, new Matrix3())

  /* step5. 判断条件 */
  if (times >= maxTime || jacobiLoopCondition(tempMatrix)) {
    return {
      eigenValues: [tempMatrix.m11, tempMatrix.m22, tempMatrix.m33], // 特征值
      eigenVectors: [
        [givensRight.m11, givensRight.m21, givensRight.m31],
        [givensRight.m12, givensRight.m22, givensRight.m32],
        [givensRight.m13, givensRight.m23, givensRight.m33]
      ] 
    }
  } else {
    return jacobi3d(tempMatrix)
  }
}

export default jacobi

计算二三维实对称矩阵的特征值、特征向量
例如

import jacobi from './alg/jacobi.js'
import { fromValuesArray } from './core/mat3.js'
const data = [
  4, 2, 2,
  2, 5, 1,
  2, 1, 6
]

const m = fromValuesArray(data)
const result = jacobi(m)

得到三个近似特征值：2.1260943732719926, 8.387618874579875, 4.486286752148127
对应三个特征向量：
[0.8292299607559941, -0.46364488168716955, -0.31210750691052425]
[0.5388204749371205, 0.5148049627855905, 0.6668195753574295]
[-0.14849298964310811, -0.7211166854099392, 0.676713054440532]

与C++、网络计算版本相差在 5% 内
算法有待优化，受限于求旋转角的精度瑕疵，争取能将特征向量逼近到1%内