Skip to content

Commit 55b1669

Browse files
tomitomi3tomitomi3
committed
コードの整理
1 parent a54a2ce commit 55b1669

File tree

3 files changed

+57
-221
lines changed

3 files changed

+57
-221
lines changed

‎LibOptimization/EasyLinearAlgebra/clsEasyMatrix.vb‎

Lines changed: 17 additions & 204 deletions
Original file line numberDiff line numberDiff line change
@@ -1419,7 +1419,7 @@
14191419
'transopose
14201420
'matP = matP.T
14211421

1422-
a.PrintValue(10)
1422+
'a.PrintValue(10)
14231423

14241424
'replace
14251425
Dim matL = New clsEasyMatrix(n, True)
@@ -1486,103 +1486,16 @@
14861486
Return ret
14871487
End Function
14881488

1489-
''' <summary>
1490-
''' To Tridiagonal matrix using Householder transform 三重対角行列
1491-
''' </summary>
1492-
''' <returns></returns>
1493-
Public Function HouseholderTransformation() As clsEasyMatrix
1494-
'ref
1495-
'ハウスホルダー変換
1496-
'http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?Numerical%20Calculation
1497-
1498-
Dim tempMat = New clsEasyMatrix(Me)
1499-
Dim n As Integer = tempMat.RowCount
1500-
Dim u = New clsEasyVector(n)
1501-
Dim v = New clsEasyVector(n)
1502-
Dim q = New clsEasyVector(n)
1503-
1504-
For k As Integer = 0 To n - 2
1505-
's
1506-
Dim s As Double = 0.0
1507-
For i = k + 1 To n - 1
1508-
s += tempMat(i)(k) * tempMat(i)(k)
1509-
Next
1510-
1511-
' | -
1512-
' |a10 -
1513-
' | a21
1514-
Dim temp As Double = tempMat(k + 1)(k)
1515-
If temp >= 0 Then
1516-
s = -Math.Sqrt(s)
1517-
Else
1518-
s = Math.Sqrt(s)
1519-
End If
1520-
1521-
' |x-y|
1522-
Dim alpha = Math.Sqrt(2.0 * s * (s - temp))
1523-
If clsMathUtil.IsCloseToZero(alpha) = True Then
1524-
Continue For
1525-
End If
1526-
1527-
'u
1528-
u(k + 1) = (temp - s) / alpha
1529-
For i = k + 2 To n - 1
1530-
u(i) = tempMat(i)(k) / alpha
1531-
Next
1532-
1533-
'Au
1534-
q(k) = alpha / 2.0
1535-
For i = k + 1 To n - 1
1536-
q(i) = 0.0
1537-
For j = k + 1 To n - 1
1538-
q(i) += tempMat(i)(j) * u(j)
1539-
Next
1540-
Next
1541-
1542-
'v=2(Au-uu^T(Au))
1543-
alpha = 0.0
1544-
'uu^T
1545-
For i = k + 1 To n - 1
1546-
alpha += u(i) * q(i)
1547-
Next
1548-
v(k) = 2.0 * q(k)
1549-
For i = k + 1 To n - 1
1550-
v(i) = 2.0 * (q(i) - alpha * u(i))
1551-
Next
1552-
1553-
'A = PAP
1554-
' | - a02
1555-
' | -
1556-
' |a20 -
1557-
tempMat(k)(k + 1) = s
1558-
tempMat(k + 1)(k) = s
1559-
For i = k + 2 To n - 1
1560-
tempMat(k)(i) = 0.0
1561-
tempMat(i)(k) = 0.0
1562-
Next
1563-
For i = k + 1 To n - 1
1564-
tempMat(i)(i) = tempMat(i)(i) - 2.0 * u(i) * v(i)
1565-
For j = i + 1 To n - 1
1566-
Dim tempVal = tempMat(i)(j) - u(i) * v(j) - v(i) * u(j)
1567-
tempMat(i)(j) = tempVal
1568-
tempMat(j)(i) = tempVal
1569-
Next
1570-
Next
1571-
Next
1572-
1573-
Return tempMat
1574-
End Function
1575-
15761489
''' <summary>
15771490
''' Eigen decomposition using Jacobi Method. for Symmetric Matrix.
15781491
''' Memo: A = V*D*V−1, D is diag(eigen value1 ... eigen valueN), V is eigen vectors. V is orthogonal matrix.
15791492
''' </summary>
15801493
''' <param name="Iteration">default:1000</param>
1581-
''' <param name="Conversion">default:1.0e-16</param>
1494+
''' <param name="Conversion">default:1.0e-15</param>
15821495
''' <param name="IsSort">descent sort by EigenValue default:true</param>
15831496
''' <returns></returns>
15841497
Public Function Eigen(Optional ByVal Iteration As Integer = 1000,
1585-
Optional ByVal Conversion As Double = 0.0000000000000001,
1498+
Optional ByVal Conversion As Double = 0.000000000000001,
15861499
Optional ByVal IsSort As Boolean = True) As Eigen
15871500
'固有値、固有ベクトルを求める方針
15881501
' 反復計算によって求める。計算しやすい行列に変換
@@ -1675,7 +1588,7 @@
16751588
'sort by Eigen value
16761589
Dim eigenValue = retEigenMat.ToDiagonalVector()
16771590
If IsSort = True Then
1678-
clsMathUtil.EigenSort(eigenValue, rotate)
1591+
clsMathUtil.EigenSort(eigenValue, rotate, True)
16791592
End If
16801593

16811594
Return New Eigen(eigenValue, rotate, isConversion)
@@ -1684,12 +1597,12 @@
16841597
''' <summary>
16851598
''' Eigen
16861599
''' </summary>
1687-
''' <param name="Iteration"></param>
1688-
''' <param name="Conversion"></param>
1689-
''' <param name="IsSort"></param>
1600+
''' <param name="Iteration">default:1000</param>
1601+
''' <param name="Conversion">default:1.0e-15</param>
1602+
''' <param name="IsSort">descent sort by EigenValue default:true</param>
16901603
''' <returns></returns>
16911604
Public Function Eigen2(Optional ByVal Iteration As Integer = 1000,
1692-
Optional ByVal Conversion As Double = 0.0000000000000001,
1605+
Optional ByVal Conversion As Double = 0.000000000000001,
16931606
Optional ByVal IsSort As Boolean = True) As Eigen
16941607
If Me.IsSquare() = False Then
16951608
Throw New clsException(clsException.Series.DifferRowNumberAndCollumnNumber)
@@ -1838,7 +1751,7 @@
18381751

18391752
'sort by Eigen value
18401753
If IsSort = True Then
1841-
clsMathUtil.EigenSort(eigenValue, eigenVector)
1754+
clsMathUtil.EigenSort(eigenValue, eigenVector, False)
18421755
End If
18431756

18441757
Return New Eigen(eigenValue, eigenVector, True)
@@ -1847,9 +1760,9 @@
18471760
''' <summary>
18481761
''' Eigen
18491762
''' </summary>
1850-
''' <param name="Iteration"></param>
1851-
''' <param name="Conversion"></param>
1852-
''' <param name="IsSort"></param>
1763+
''' <param name="Iteration">default:1000</param>
1764+
''' <param name="Conversion">default:1.0e-15</param>
1765+
''' <param name="IsSort">descent sort by EigenValue default:true</param>
18531766
''' <returns></returns>
18541767
Public Function Eigen3(Optional ByVal Iteration As Integer = 1000,
18551768
Optional ByVal Conversion As Double = 0.000000000000001,
@@ -1869,7 +1782,6 @@
18691782
While (m > 1)
18701783
Dim dVal = a(m)(m - 1)
18711784
If Math.Abs(dVal) < Conversion Then
1872-
a.PrintValue()
18731785
m -= 1
18741786
End If
18751787

@@ -1989,7 +1901,7 @@
19891901

19901902
'sort by Eigen value
19911903
If IsSort = True Then
1992-
clsMathUtil.EigenSort(eigenValues, eigenVectors)
1904+
clsMathUtil.EigenSort(eigenValues, eigenVectors, False)
19931905
End If
19941906

19951907
Return New Eigen(eigenValues, eigenVectors, True)
@@ -1998,9 +1910,9 @@
19981910
''' <summary>
19991911
''' Eigen
20001912
''' </summary>
2001-
''' <param name="Iteration"></param>
2002-
''' <param name="Conversion"></param>
2003-
''' <param name="IsSort"></param>
1913+
''' <param name="Iteration">default:1000</param>
1914+
''' <param name="Conversion">default:1.0e-15</param>
1915+
''' <param name="IsSort">descent sort by EigenValue default:true</param>
20041916
''' <returns></returns>
20051917
Public Function Eigen4(Optional ByVal Iteration As Integer = 1000,
20061918
Optional ByVal Conversion As Double = 0.000000000000001,
@@ -2020,7 +1932,6 @@
20201932
While (m > 1)
20211933
Dim dVal = a(m)(m - 1)
20221934
If Math.Abs(dVal) < Conversion Then
2023-
a.PrintValue()
20241935
m -= 1
20251936
End If
20261937

@@ -2133,7 +2044,7 @@
21332044
Next
21342045
v(k) = temp / ludecomp(k)(k)
21352046
Next
2136-
v.PrintValue()
2047+
'v.PrintValue()
21372048

21382049
mu = v.InnerProduct(y)
21392050
v2 = v.NormL2()
@@ -2287,104 +2198,6 @@
22872198
Next
22882199
Return a
22892200
End Function
2290-
2291-
''' <summary>
2292-
''' Householder Transformation
2293-
''' 非対称行列をハウスホルダー変換 → ヘッセンベルグ行列
2294-
''' 対称行列をハウスホルダー変換 → 三重対角行列
2295-
''' </summary>
2296-
''' <returns></returns>
2297-
Public Function HouseholderTransformationForQR() As clsEasyMatrix
2298-
'ref
2299-
'ハウスホルダー変換
2300-
'http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?Numerical%20Calculation
2301-
2302-
Dim a = New clsEasyMatrix(Me)
2303-
Dim n = a.Count
2304-
2305-
Dim b = New clsEasyMatrix(n)
2306-
Dim p = New clsEasyMatrix(n)
2307-
Dim q = New clsEasyMatrix(n)
2308-
Dim u = New clsEasyVector(n)
2309-
2310-
For k As Integer = 0 To n - 3
2311-
's
2312-
Dim s As Double = 0.0
2313-
For i = k + 1 To n - 1
2314-
s += a(i)(k) * a(i)(k)
2315-
Next
2316-
Dim tempVal = a(k + 1)(k)
2317-
If tempVal >= 0 Then
2318-
s = -Math.Sqrt(s)
2319-
Else
2320-
s = Math.Sqrt(s)
2321-
End If
2322-
2323-
' |x-y|
2324-
Dim alpha = Math.Sqrt(2.0 * s * (s - tempVal))
2325-
If clsMathUtil.IsCloseToZero(alpha) = True Then
2326-
Continue For
2327-
End If
2328-
2329-
'u
2330-
u(k + 1) = (tempVal - s) / alpha
2331-
For i = k + 2 To n - 1
2332-
u(i) = a(i)(k) / alpha
2333-
Next
2334-
2335-
'P
2336-
For i = k + 1 To n - 1
2337-
For j = i To n - 1
2338-
If j = i Then
2339-
p(i)(j) = 1.0 - 2.0 * u(i) * u(i)
2340-
Else
2341-
p(i)(j) = -2.0 * u(i) * u(j)
2342-
p(j)(i) = p(i)(j)
2343-
End If
2344-
Next
2345-
Next
2346-
2347-
'PA
2348-
For i = k + 1 To n - 1
2349-
For j = k + 1 To n - 1
2350-
q(i)(j) = 0.0
2351-
For m = k + 1 To n - 1
2352-
q(i)(j) += p(i)(m) * a(m)(j)
2353-
Next
2354-
Next
2355-
Next
2356-
2357-
'A = PAP^T
2358-
For i = 0 To k
2359-
b(i)(k) = a(i)(k)
2360-
Next
2361-
b(k + 1)(k) = s
2362-
For i = k + 2 To n - 1
2363-
b(i)(k) = 0.0
2364-
Next
2365-
For j = k + 1 To n - 1
2366-
For i = 0 To k
2367-
b(i)(j) = 0.0
2368-
For m = k + 1 To n - 1
2369-
b(i)(j) += a(i)(m) * p(j)(m)
2370-
Next
2371-
Next
2372-
For i = k + 1 To n - 1
2373-
b(i)(j) = 0.0
2374-
For m = k + 1 To n - 1
2375-
b(i)(j) += q(i)(m) * p(j)(m)
2376-
Next
2377-
Next
2378-
Next
2379-
2380-
For i = 0 To n - 1
2381-
For j = 0 To n - 1
2382-
a(i)(j) = b(i)(j)
2383-
Next
2384-
Next
2385-
Next
2386-
Return a
2387-
End Function
23882201
#End Region
23892202

23902203
#Region "Public Shared"

‎LibOptimization/EasyLinearAlgebra/clsMathUtil.vb‎

Lines changed: 25 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -244,26 +244,40 @@
244244
''' </summary>
245245
''' <param name="eigenValue"></param>
246246
''' <param name="eigenVector"></param>
247-
Public Shared Sub EigenSort(ByRef eigenValue As clsEasyVector, ByRef eigenVector As clsEasyMatrix)
247+
''' <param name="isColOrder"></param>
248+
Public Shared Sub EigenSort(ByRef eigenValue As clsEasyVector, ByRef eigenVector As clsEasyMatrix, ByVal isColOrder As Boolean)
248249
Dim n = eigenValue.Count
249250
Dim colSwapInfo = New List(Of ValueDescSort)
250251
For i As Integer = 0 To n - 1
251252
colSwapInfo.Add(New ValueDescSort(eigenValue(i), i))
252253
Next
253254
colSwapInfo.Sort()
254255

255-
Dim newEigenVector = New clsEasyMatrix(n)
256-
For j As Integer = 0 To n - 1
257-
'eigen value
258-
eigenValue(j) = colSwapInfo(j).v
256+
If isColOrder = True Then
257+
Dim newEigenVector = New clsEasyMatrix(n)
258+
For j As Integer = 0 To n - 1
259+
'eigen value
260+
eigenValue(j) = colSwapInfo(j).v
259261

260-
'eigen vector
261-
Dim k = colSwapInfo(j).idx
262-
For i As Integer = 0 To n - 1
263-
newEigenVector(i)(j) = eigenVector(i)(k)
262+
'eigen vector
263+
Dim k = colSwapInfo(j).idx
264+
For i As Integer = 0 To n - 1
265+
newEigenVector(i)(j) = eigenVector(i)(k)
266+
Next
264267
Next
265-
Next
266-
eigenVector = newEigenVector
268+
eigenVector = newEigenVector
269+
Else
270+
Dim newEigenVector = New clsEasyMatrix(n)
271+
For i = 0 To n - 1
272+
'eigen value
273+
eigenValue(i) = colSwapInfo(i).v
274+
275+
'eigen vector
276+
Dim k = colSwapInfo(i).idx
277+
newEigenVector(i) = eigenVector(k)
278+
Next
279+
eigenVector = newEigenVector
280+
End If
267281
End Sub
268282

269283
''' <summary>

0 commit comments

Comments
 (0)