| 478 |
478 |
end
|
| 479 |
479 |
|
| 480 |
480 |
#
|
| 481 |
|
# Returns +true+ if this is a regular matrix.
|
|
481 |
# Returns +true+ if this is a regular (i.e. non-singular) matrix.
|
| 482 |
482 |
#
|
| 483 |
483 |
def regular?
|
| 484 |
|
square? and rank == column_size
|
|
484 |
not singular?
|
| 485 |
485 |
end
|
| 486 |
486 |
|
| 487 |
487 |
#
|
| 488 |
|
# Returns +true+ is this is a singular (i.e. non-regular) matrix.
|
|
488 |
# Returns +true+ is this is a singular matrix.
|
| 489 |
489 |
#
|
| 490 |
490 |
def singular?
|
| 491 |
|
not regular?
|
|
491 |
determinant == 0
|
| 492 |
492 |
end
|
| 493 |
493 |
|
| 494 |
494 |
#
|
| ... | ... | |
| 740 |
740 |
|
| 741 |
741 |
#
|
| 742 |
742 |
# Returns the determinant of the matrix.
|
| 743 |
|
# This method's algorithm is Gaussian elimination method
|
| 744 |
|
# and using Numeric#quo(). Beware that using Float values, with their
|
| 745 |
|
# usual lack of precision, can affect the value returned by this method. Use
|
| 746 |
|
# Rational values or Matrix#det_e instead if this is important to you.
|
|
743 |
#
|
|
744 |
# Beware that using Float values can yield erroneous results
|
|
745 |
# because of their lack of precision.
|
|
746 |
# Consider using exact types like Rational or BigDecimal instead.
|
| 747 |
747 |
#
|
| 748 |
748 |
# Matrix[[7,6], [3,9]].determinant
|
| 749 |
|
# => 45.0
|
|
749 |
# => 45
|
| 750 |
750 |
#
|
| 751 |
751 |
def determinant
|
| 752 |
752 |
Matrix.Raise ErrDimensionMismatch unless square?
|
|
753 |
m = @rows
|
|
754 |
case row_size
|
|
755 |
# Up to 4x4, give result using Laplacian expansion by minors.
|
|
756 |
# This will typically be faster, as well as giving good results
|
|
757 |
# in case of Floats
|
|
758 |
when 0
|
|
759 |
+1
|
|
760 |
when 1
|
|
761 |
+ m[0][0]
|
|
762 |
when 2
|
|
763 |
+ m[0][0] * m[1][1] - m[0][1] * m[1][0]
|
|
764 |
when 3
|
|
765 |
m0 = m[0]; m1 = m[1]; m2 = m[2]
|
|
766 |
+ m0[0] * m1[1] * m2[2] - m0[0] * m1[2] * m2[1] \
|
|
767 |
- m0[1] * m1[0] * m2[2] + m0[1] * m1[2] * m2[0] \
|
|
768 |
+ m0[2] * m1[0] * m2[1] - m0[2] * m1[1] * m2[0]
|
|
769 |
when 4
|
|
770 |
m0 = m[0]; m1 = m[1]; m2 = m[2]; m3 = m[3]
|
|
771 |
+ m0[0] * m1[1] * m2[2] * m3[3] - m0[0] * m1[1] * m2[3] * m3[2] \
|
|
772 |
- m0[0] * m1[2] * m2[1] * m3[3] + m0[0] * m1[2] * m2[3] * m3[1] \
|
|
773 |
+ m0[0] * m1[3] * m2[1] * m3[2] - m0[0] * m1[3] * m2[2] * m3[1] \
|
|
774 |
- m0[1] * m1[0] * m2[2] * m3[3] + m0[1] * m1[0] * m2[3] * m3[2] \
|
|
775 |
+ m0[1] * m1[2] * m2[0] * m3[3] - m0[1] * m1[2] * m2[3] * m3[0] \
|
|
776 |
- m0[1] * m1[3] * m2[0] * m3[2] + m0[1] * m1[3] * m2[2] * m3[0] \
|
|
777 |
+ m0[2] * m1[0] * m2[1] * m3[3] - m0[2] * m1[0] * m2[3] * m3[1] \
|
|
778 |
- m0[2] * m1[1] * m2[0] * m3[3] + m0[2] * m1[1] * m2[3] * m3[0] \
|
|
779 |
+ m0[2] * m1[3] * m2[0] * m3[1] - m0[2] * m1[3] * m2[1] * m3[0] \
|
|
780 |
- m0[3] * m1[0] * m2[1] * m3[2] + m0[3] * m1[0] * m2[2] * m3[1] \
|
|
781 |
+ m0[3] * m1[1] * m2[0] * m3[2] - m0[3] * m1[1] * m2[2] * m3[0] \
|
|
782 |
- m0[3] * m1[2] * m2[0] * m3[1] + m0[3] * m1[2] * m2[1] * m3[0]
|
|
783 |
else
|
|
784 |
# For bigger matrices, use an efficient and general algorithm.
|
|
785 |
# Currently, we use the Gauss-Bareiss algorithm
|
|
786 |
determinant_bareiss
|
|
787 |
end
|
|
788 |
end
|
|
789 |
alias_method :det, :determinant
|
| 753 |
790 |
|
|
791 |
#
|
|
792 |
# Private. Use Matrix#determinant
|
|
793 |
#
|
|
794 |
# Returns the determinant of the matrix, using
|
|
795 |
# Bareiss' multistep integer-preserving gaussian elimination.
|
|
796 |
# It has the same computational cost order O(n^3) as standard Gaussian elimination.
|
|
797 |
# Intermediate results are fraction free and of lower complexity.
|
|
798 |
# A matrix of Integers will have thus intermediate results that are also Integers,
|
|
799 |
# with smaller bignums (if any), while a matrix of Float will usually have more
|
|
800 |
# precise intermediate results.
|
|
801 |
#
|
|
802 |
def determinant_bareiss
|
| 754 |
803 |
size = row_size
|
|
804 |
last = size - 1
|
| 755 |
805 |
a = to_a
|
| 756 |
|
|
| 757 |
|
det = 1
|
|
806 |
no_pivot = Proc.new{ return 0 }
|
|
807 |
sign = +1
|
|
808 |
pivot = 1
|
| 758 |
809 |
size.times do |k|
|
| 759 |
|
if (akk = a[k][k]) == 0
|
| 760 |
|
i = (k+1 ... size).find {|ii|
|
| 761 |
|
a[ii][k] != 0
|
|
810 |
previous_pivot = pivot
|
|
811 |
if (pivot = a[k][k]) == 0
|
|
812 |
switch = (k+1 ... size).find(no_pivot) {|row|
|
|
813 |
a[row][k] != 0
|
| 762 |
814 |
}
|
| 763 |
|
return 0 if i.nil?
|
| 764 |
|
a[i], a[k] = a[k], a[i]
|
| 765 |
|
akk = a[k][k]
|
| 766 |
|
det *= -1
|
|
815 |
a[switch], a[k] = a[k], a[switch]
|
|
816 |
pivot = a[k][k]
|
|
817 |
sign = -sign
|
| 767 |
818 |
end
|
| 768 |
|
|
| 769 |
|
(k + 1 ... size).each do |ii|
|
| 770 |
|
q = a[ii][k].quo(akk)
|
| 771 |
|
(k + 1 ... size).each do |j|
|
| 772 |
|
a[ii][j] -= a[k][j] * q
|
|
819 |
(k+1).upto(last) do |i|
|
|
820 |
ai = a[i]
|
|
821 |
(k+1).upto(last) do |j|
|
|
822 |
ai[j] = (pivot * ai[j] - ai[k] * a[k][j]) / previous_pivot
|
| 773 |
823 |
end
|
| 774 |
824 |
end
|
| 775 |
|
det *= akk
|
| 776 |
825 |
end
|
| 777 |
|
det
|
|
826 |
sign * pivot
|
| 778 |
827 |
end
|
| 779 |
|
alias det determinant
|
|
828 |
private :determinant_bareiss
|
| 780 |
829 |
|
| 781 |
830 |
#
|
| 782 |
|
# Returns the determinant of the matrix.
|
| 783 |
|
# This method's algorithm is Gaussian elimination method.
|
| 784 |
|
# This method uses Euclidean algorithm. If all elements are integer,
|
| 785 |
|
# really exact value. But, if an element is a float, can't return
|
| 786 |
|
# exact value.
|
| 787 |
|
#
|
| 788 |
|
# Matrix[[7,6], [3,9]].determinant
|
| 789 |
|
# => 63
|
|
831 |
# deprecated; use Matrix#determinant
|
| 790 |
832 |
#
|
| 791 |
833 |
def determinant_e
|
| 792 |
|
Matrix.Raise ErrDimensionMismatch unless square?
|
| 793 |
|
|
| 794 |
|
size = row_size
|
| 795 |
|
a = to_a
|
| 796 |
|
|
| 797 |
|
det = 1
|
| 798 |
|
size.times do |k|
|
| 799 |
|
if a[k][k].zero?
|
| 800 |
|
i = (k+1 ... size).find {|ii|
|
| 801 |
|
a[ii][k] != 0
|
| 802 |
|
}
|
| 803 |
|
return 0 if i.nil?
|
| 804 |
|
a[i], a[k] = a[k], a[i]
|
| 805 |
|
det *= -1
|
| 806 |
|
end
|
| 807 |
|
|
| 808 |
|
(k + 1 ... size).each do |ii|
|
| 809 |
|
q = a[ii][k].quo(a[k][k])
|
| 810 |
|
(k ... size).each do |j|
|
| 811 |
|
a[ii][j] -= a[k][j] * q
|
| 812 |
|
end
|
| 813 |
|
unless a[ii][k].zero?
|
| 814 |
|
a[ii], a[k] = a[k], a[ii]
|
| 815 |
|
det *= -1
|
| 816 |
|
redo
|
| 817 |
|
end
|
| 818 |
|
end
|
| 819 |
|
det *= a[k][k]
|
| 820 |
|
end
|
| 821 |
|
det
|
|
834 |
warn "#{caller(1)[0]}: warning: Matrix#determinant_e is deprecated; use #determinant"
|
|
835 |
rank
|
| 822 |
836 |
end
|
| 823 |
837 |
alias det_e determinant_e
|
| 824 |
838 |
|
| 825 |
839 |
#
|
| 826 |
|
# Returns the rank of the matrix. Beware that using Float values,
|
| 827 |
|
# probably return faild value. Use Rational values or Matrix#rank_e
|
| 828 |
|
# for getting exact result.
|
|
840 |
# Returns the rank of the matrix.
|
|
841 |
# Beware that using Float values can yield erroneous results
|
|
842 |
# because of their lack of precision.
|
|
843 |
# Consider using exact types like Rational or BigDecimal instead.
|
| 829 |
844 |
#
|
| 830 |
845 |
# Matrix[[7,6], [3,9]].rank
|
| 831 |
846 |
# => 2
|
| 832 |
847 |
#
|
| 833 |
848 |
def rank
|
| 834 |
|
if column_size > row_size
|
| 835 |
|
a = transpose.to_a
|
| 836 |
|
a_column_size = row_size
|
| 837 |
|
a_row_size = column_size
|
| 838 |
|
else
|
| 839 |
|
a = to_a
|
| 840 |
|
a_column_size = column_size
|
| 841 |
|
a_row_size = row_size
|
| 842 |
|
end
|
|
849 |
# We currently use Bareiss' multistep integer-preserving gaussian elimination
|
|
850 |
# (see comments on determinant)
|
|
851 |
a = to_a
|
|
852 |
last_column = column_size - 1
|
|
853 |
last_row = row_size - 1
|
| 843 |
854 |
rank = 0
|
| 844 |
|
a_column_size.times do |k|
|
| 845 |
|
if (akk = a[k][k]) == 0
|
| 846 |
|
i = (k+1 ... a_row_size).find {|ii|
|
| 847 |
|
a[ii][k] != 0
|
| 848 |
|
}
|
| 849 |
|
if i
|
| 850 |
|
a[i], a[k] = a[k], a[i]
|
| 851 |
|
akk = a[k][k]
|
| 852 |
|
else
|
| 853 |
|
i = (k+1 ... a_column_size).find {|ii|
|
| 854 |
|
a[k][ii] != 0
|
| 855 |
|
}
|
| 856 |
|
next if i.nil?
|
| 857 |
|
(k ... a_column_size).each do |j|
|
| 858 |
|
a[j][k], a[j][i] = a[j][i], a[j][k]
|
| 859 |
|
end
|
| 860 |
|
akk = a[k][k]
|
| 861 |
|
end
|
|
855 |
pivot_row = 0
|
|
856 |
previous_pivot = 1
|
|
857 |
0.upto(last_column) do |k|
|
|
858 |
switch_row = (pivot_row .. last_row).find {|row|
|
|
859 |
a[row][k] != 0
|
|
860 |
}
|
|
861 |
if switch_row
|
|
862 |
a[switch_row], a[pivot_row] = a[pivot_row], a[switch_row] unless pivot_row == switch_row
|
|
863 |
pivot = a[pivot_row][k]
|
|
864 |
(pivot_row+1).upto(last_row) do |i|
|
|
865 |
ai = a[i]
|
|
866 |
(k+1).upto(last_column) do |j|
|
|
867 |
ai[j] = (pivot * ai[j] - ai[k] * a[pivot_row][j]) / previous_pivot
|
|
868 |
end
|
|
869 |
end
|
|
870 |
pivot_row += 1
|
|
871 |
previous_pivot = pivot
|
| 862 |
872 |
end
|
| 863 |
|
|
| 864 |
|
(k + 1 ... a_row_size).each do |ii|
|
| 865 |
|
q = a[ii][k].quo(akk)
|
| 866 |
|
(k + 1... a_column_size).each do |j|
|
| 867 |
|
a[ii][j] -= a[k][j] * q
|
| 868 |
|
end
|
| 869 |
|
end
|
| 870 |
|
rank += 1
|
| 871 |
873 |
end
|
| 872 |
|
return rank
|
|
874 |
pivot_row
|
| 873 |
875 |
end
|
| 874 |
876 |
|
| 875 |
877 |
#
|
| 876 |
|
# Returns the rank of the matrix. This method uses Euclidean
|
| 877 |
|
# algorithm. If all elements are integer, really exact value. But, if
|
| 878 |
|
# an element is a float, can't return exact value.
|
| 879 |
|
#
|
| 880 |
|
# Matrix[[7,6], [3,9]].rank
|
| 881 |
|
# => 2
|
|
878 |
# deprecated; use Matrix#rank
|
| 882 |
879 |
#
|
| 883 |
880 |
def rank_e
|
| 884 |
|
a = to_a
|
| 885 |
|
a_column_size = column_size
|
| 886 |
|
a_row_size = row_size
|
| 887 |
|
pi = 0
|
| 888 |
|
a_column_size.times do |j|
|
| 889 |
|
if i = (pi ... a_row_size).find{|i0| !a[i0][j].zero?}
|
| 890 |
|
if i != pi
|
| 891 |
|
a[pi], a[i] = a[i], a[pi]
|
| 892 |
|
end
|
| 893 |
|
(pi + 1 ... a_row_size).each do |k|
|
| 894 |
|
q = a[k][j].quo(a[pi][j])
|
| 895 |
|
(pi ... a_column_size).each do |j0|
|
| 896 |
|
a[k][j0] -= q * a[pi][j0]
|
| 897 |
|
end
|
| 898 |
|
if k > pi && !a[k][j].zero?
|
| 899 |
|
a[k], a[pi] = a[pi], a[k]
|
| 900 |
|
redo
|
| 901 |
|
end
|
| 902 |
|
end
|
| 903 |
|
pi += 1
|
| 904 |
|
end
|
| 905 |
|
end
|
| 906 |
|
pi
|
|
881 |
warn "#{caller(1)[0]}: warning: Matrix#rank_e is deprecated; use #rank"
|
|
882 |
rank
|
| 907 |
883 |
end
|
| 908 |
884 |
|
| 909 |
885 |
|