/
arithmetic.jl
132 lines (94 loc) · 3.33 KB
/
arithmetic.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
import Base: +, -, *, /, \
# =========================
# Addition operations
# =========================
+(M1::IntervalMatrix, M2::IntervalMatrix) = IntervalMatrix(M1.mat + M2.mat)
-(M1::IntervalMatrix, M2::IntervalMatrix) = IntervalMatrix(M1.mat - M2.mat)
+(x::Interval, M::IntervalMatrix) = IntervalMatrix(x .+ M.mat)
+(M::IntervalMatrix, x::Interval) = IntervalMatrix(x .+ M.mat)
+(x::Number, M::IntervalMatrix) = interval(x) + M
+(M::IntervalMatrix, x::Number) = interval(x) + M
+(M1::IntervalMatrix, M2::AbstractMatrix) = IntervalMatrix(M1.mat + M2)
+(M1::AbstractMatrix, M2::IntervalMatrix) = IntervalMatrix(M1 + M2.mat)
-(M1::IntervalMatrix, M2::AbstractMatrix) = IntervalMatrix(M1.mat - M2)
-(M1::AbstractMatrix, M2::IntervalMatrix) = IntervalMatrix(M1 - M2.mat)
# =========================
# Multiplication operations
# =========================
# matrix-matrix multiplication is defined in a separate file
*(x::Interval, M::IntervalMatrix) = IntervalMatrix(x .* M.mat)
*(M::IntervalMatrix, x::Interval) = IntervalMatrix(x .* M.mat)
*(x::Number, M::IntervalMatrix) = interval(x) * M
*(M::IntervalMatrix, x::Number) = interval(x) * M
/(M::IntervalMatrix, x::Number) = IntervalMatrix(M ./ x)
# left-division methods to avoid a stack overflow with the default behavior
# (there exist more precise approaches but are currently not implemented here)
\(M1::IntervalMatrix, M2::IntervalMatrix) = IntervalMatrix(M1.mat \ M2.mat)
\(M1::IntervalMatrix, M2::AbstractMatrix) = IntervalMatrix(M1.mat \ M2)
\(M1::AbstractMatrix, M2::IntervalMatrix) = IntervalMatrix(M1 \ M2.mat)
"""
square(A::IntervalMatrix)
Compute the square of an interval matrix.
### Input
- `A` -- interval matrix
### Output
An interval matrix equivalent to `A * A`.
### Algorithm
We follow [1, Section 6].
[1] Kosheleva, Kreinovich, Mayer, Nguyen. Computing the cube of an interval
matrix is NP-hard. SAC 2005.
"""
function square(A::IntervalMatrix)
B = similar(A)
n = checksquare(A)
# case i == j
@inbounds for j in 1:n
B[j, j] = pow(A[j, j], 2)
for k in 1:n
k == j && continue
B[j, j] += A[j, k] * A[k, j]
end
end
# case i ≠ j
@inbounds for j in 1:n
for i in 1:n
i == j && continue
B[i, j] = A[i, j] * (A[j, j] + A[i, i])
for k in 1:n
(k == i || k == j) && continue
B[i, j] += A[i, k] * A[k, j]
end
end
end
return B
end
"""
scale(A::IntervalMatrix{T}, α::T) where {T}
Return a new interval matrix whose entries are scaled by the given factor.
### Input
- `A` -- interval matrix
- `α` -- scaling factor
### Output
A new matrix `B` of the same shape as `A` such that `B[i, j] = α*A[i, j]` for
each `i` and `j`.
### Notes
See `scale!` for the in-place version of this function.
"""
function scale(A::IntervalMatrix{T}, α::T) where {T}
return scale!(copy(A), α)
end
"""
scale!(A::IntervalMatrix{T}, α::T) where {T}
Modifies the given interval matrix, scaling its entries by the given factor.
### Input
- `A` -- interval matrix
- `α` -- scaling factor
### Output
The matrix `A` such that for each `i` and `j`, the new value of `A[i, j]` is
`α*A[i, j]`.
### Notes
This is the in-place version of `scale`.
"""
function scale!(A::IntervalMatrix{T}, α::T) where {T}
return map!(x -> α * x, A, A)
end