class Matrix::LUP分解

对于一个 m×n 的矩阵 A，其中 m >= n，LU 分解是一个 m×n 的单位下三角矩阵 L，一个 n×n 的上三角矩阵 U，以及一个 m×m 的置换矩阵 P，使得 L*U = P*A。如果 m < n，则 L 是 m×m，U 是 m×n。

带旋转的 LUP 分解总是存在的，即使矩阵是奇异的，因此构造函数永远不会失败。LU 分解的主要用途是求解线性方程组。如果 singular? 返回 true，这将失败。

属性

pivots[R]

返回旋转索引

公共类方法

new(a) 点击以切换源代码

# File matrix-0.4.2/lib/matrix/lup_decomposition.rb, line 154
def initialize a
  raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix)
  # Use a "left-looking", dot-product, Crout/Doolittle algorithm.
  @lu = a.to_a
  @row_count = a.row_count
  @column_count = a.column_count
  @pivots = Array.new(@row_count)
  @row_count.times do |i|
     @pivots[i] = i
  end
  @pivot_sign = 1
  lu_col_j = Array.new(@row_count)

  # Outer loop.

  @column_count.times do |j|

    # Make a copy of the j-th column to localize references.

    @row_count.times do |i|
      lu_col_j[i] = @lu[i][j]
    end

    # Apply previous transformations.

    @row_count.times do |i|
      lu_row_i = @lu[i]

      # Most of the time is spent in the following dot product.

      kmax = [i, j].min
      s = 0
      kmax.times do |k|
        s += lu_row_i[k]*lu_col_j[k]
      end

      lu_row_i[j] = lu_col_j[i] -= s
    end

    # Find pivot and exchange if necessary.

    p = j
    (j+1).upto(@row_count-1) do |i|
      if (lu_col_j[i].abs > lu_col_j[p].abs)
        p = i
      end
    end
    if (p != j)
      @column_count.times do |k|
        t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t
      end
      k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k
      @pivot_sign = -@pivot_sign
    end

    # Compute multipliers.

    if (j < @row_count && @lu[j][j] != 0)
      (j+1).upto(@row_count-1) do |i|
        @lu[i][j] = @lu[i][j].quo(@lu[j][j])
      end
    end
  end
end

公共实例方法

det() 点击以切换源代码

返回 A 的行列式，从因式分解中高效计算。

# File matrix-0.4.2/lib/matrix/lup_decomposition.rb, line 79
def det
  if (@row_count != @column_count)
    raise Matrix::ErrDimensionMismatch
  end
  d = @pivot_sign
  @column_count.times do |j|
    d *= @lu[j][j]
  end
  d
end

别名为: determinant

determinant()

别名为: det

l() 点击以切换源代码

# File matrix-0.4.2/lib/matrix/lup_decomposition.rb, line 22
def l
  Matrix.build(@row_count, [@column_count, @row_count].min) do |i, j|
    if (i > j)
      @lu[i][j]
    elsif (i == j)
      1
    else
      0
    end
  end
end

p() 点击以切换源代码

返回置换矩阵 P

# File matrix-0.4.2/lib/matrix/lup_decomposition.rb, line 48
def p
  rows = Array.new(@row_count){Array.new(@row_count, 0)}
  @pivots.each_with_index{|p, i| rows[i][p] = 1}
  Matrix.send :new, rows, @row_count
end

singular?() 点击以切换源代码

如果 U，以及 A，是奇异的，则返回 true。

# File matrix-0.4.2/lib/matrix/lup_decomposition.rb, line 67
def singular?
  @column_count.times do |j|
    if (@lu[j][j] == 0)
      return true
    end
  end
  false
end

solve(b) 点击以切换源代码

返回 m，使得 A*m = b，或者等价地，使得 L*U*m = P*b b 可以是一个 Matrix 或者一个 Vector

# File matrix-0.4.2/lib/matrix/lup_decomposition.rb, line 95
def solve b
  if (singular?)
    raise Matrix::ErrNotRegular, "Matrix is singular."
  end
  if b.is_a? Matrix
    if (b.row_count != @row_count)
      raise Matrix::ErrDimensionMismatch
    end

    # Copy right hand side with pivoting
    nx = b.column_count
    m = @pivots.map{|row| b.row(row).to_a}

    # Solve L*Y = P*b
    @column_count.times do |k|
      (k+1).upto(@column_count-1) do |i|
        nx.times do |j|
          m[i][j] -= m[k][j]*@lu[i][k]
        end
      end
    end
    # Solve U*m = Y
    (@column_count-1).downto(0) do |k|
      nx.times do |j|
        m[k][j] = m[k][j].quo(@lu[k][k])
      end
      k.times do |i|
        nx.times do |j|
          m[i][j] -= m[k][j]*@lu[i][k]
        end
      end
    end
    Matrix.send :new, m, nx
  else # same algorithm, specialized for simpler case of a vector
    b = convert_to_array(b)
    if (b.size != @row_count)
      raise Matrix::ErrDimensionMismatch
    end

    # Copy right hand side with pivoting
    m = b.values_at(*@pivots)

    # Solve L*Y = P*b
    @column_count.times do |k|
      (k+1).upto(@column_count-1) do |i|
        m[i] -= m[k]*@lu[i][k]
      end
    end
    # Solve U*m = Y
    (@column_count-1).downto(0) do |k|
      m[k] = m[k].quo(@lu[k][k])
      k.times do |i|
        m[i] -= m[k]*@lu[i][k]
      end
    end
    Vector.elements(m, false)
  end
end

to_a()

别名为: to_ary

to_ary() 点击以切换源代码

返回数组中的 L、U、P

# File matrix-0.4.2/lib/matrix/lup_decomposition.rb, line 56
def to_ary
  [l, u, p]
end

别名为: to_a

u() 点击以切换源代码

返回上三角因子 U

# File matrix-0.4.2/lib/matrix/lup_decomposition.rb, line 36
def u
  Matrix.build([@column_count, @row_count].min, @column_count) do |i, j|
    if (i <= j)
      @lu[i][j]
    else
      0
    end
  end
end