Matrix::LUPDecomposition

class Matrix::LUPDecomposition

Parent:Object

对于m≥n的m×n矩阵A，LU分解是m×n单位下三角矩阵L，n×n上三角矩阵U和m×m置换矩阵P使得L * U = P * A。如果m <n，那么L为m乘m，U为m乘n。

即使矩阵是单数的，使用pivoting的LUP分解总是存在的，所以构造函数永远不会失败。LU分解的主要用途是求解联立线性方程组的平方系统。如果是单数，则会失败吗？返回true。

属性

pivotsR

返回旋转索引

公共类方法

new(a) Show source

# File lib/matrix/lup_decomposition.rb, line 153
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() Show source

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

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

也被命名为：行列式

determinant()

别名为：det

l() Show source

# File lib/matrix/lup_decomposition.rb, line 21
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() Show source

返回置换矩阵 P

# File lib/matrix/lup_decomposition.rb, line 47
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?() Show source

如果U和A是单数，则返回true。

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

solve(b) Show source

返回m使得A * m = b，或等价地使得L * U * m = P * b b可以是矩阵或向量

# File lib/matrix/lup_decomposition.rb, line 94
def solve b
  if (singular?)
    Matrix.Raise Matrix::ErrNotRegular, "Matrix is singular."
  end
  if b.is_a? Matrix
    if (b.row_count != @row_count)
      Matrix.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)
      Matrix.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() Show source

返回数组中的L，U，P.

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

另外别名为：to_a

u() Show source

返回上三角因子 U

# File lib/matrix/lup_decomposition.rb, line 35
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