# Efficiency of matrix data operations

## Background:

I have the following matrix is raw data,for example:

a=\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ &...& \end{matrix} \right] , b=\left[ \begin{matrix} 1 & 2 & 3 & 5 \\ 4 & 5 & 8 & 6 \\ 7 & 0 & 9 & 8 \\ 2 & 3 & 1 & 6 \\ 4 & 5 & 6 & 9 \\ 8 & 9 & 7 & 1 \\ &...&... \end{matrix} \right]

It can be discovered:(Also what I want to get through the program)

• The first line of a : (\left[ \begin{matrix} 1 & 2 & 3 \end{matrix} \right]) is a subset of the first (\left[ \begin{matrix} 1 & 2 & 3 &5 \end{matrix} \right]) and fourth (\left[ \begin{matrix} 2 &3 &1&6 \end{matrix} \right]) lines of b.
• The first line of b differs from the first line of a by “5”,
the fourth line of b differs from the first line of a by “6”.
• According to the above ideas, the result should be like:
c=\left[ \begin{matrix} 1 & 2 & 3 &1&4&5&6\\ 4 & 5 & 6 &2&5&8&9\\ 7 & 8 & 9&3&6&0&1\\ & & &...& & \end{matrix} \right]

## My code:

c = hcat(a,zeros(Int64,size(a,1),4))
for i = 1:size(a,1)
dex = Vector{Int64}()
for j = 1:size(b,1)
if a[i,:] ⊆ b[j,:]
dex = append!(dex,j)
end
end
c[i,4] = dex[1]
c[i,5] = dex[2]
end
# Determine the number of lines a belongs to b,
# here I am sure that there are only two lines in b containing a

for i = 1:size(a,1)
k1 = setdiff(b[a[i,4],:], a[i,:])
k2 = setdiff(b[a[i,5],:], a[i,:])
c[i,6] = k1[1]
c[i,7] = k2[1]
end
# Write the element different from a in b to c


## Question:

When my matrix size is large, the execution efficiency of the above code is really too low.
I want to ask if you have any good suggestions or other algorithms to achieve the same purpose.(That is to get the matrix c)
If you still need other information, please let me know in the comment, I will add it in time, thank you all.

It is not clear to me what you are trying to do.

If you want to find repeated rows in a matrix or matrices, use an O(1) lookup mechanism, eg put the rows in a Dict.

I am actually surprise this does not throw error?
append!takes two collections,dexis a collection,j is not?

edit: huh, it works on v1.0.3, I guess I am surprise I didn’t know that.
Btw, I don’t think append! need assignment

I think this is wrong, it’s really

        k1 = setdiff(b[c[i,4],:], a[i,:])
k2 = setdiff(b[c[i,5],:], a[i,:])


I really suggest that you read the Performance Tips to improve your code.

I did not try to understand what you are trying to accomplish, but just by annotating the loops to avoid bound checking and using views, I could decrease the computational burden a lot:

function test1(a,b)
c = hcat(a,zeros(Int64,size(a,1),4))
for i = 1:size(a,1)
dex = Vector{Int64}()
for j = 1:size(b,1)
if a[i,:] ⊆ b[j,:]
dex = append!(dex,j)
end
end
c[i,4] = dex[1]
c[i,5] = dex[2]
end
# Determine the number of lines a belongs to b,
# here I am sure that there are only two lines in b containing a

for i = 1:size(a,1)
k1 = setdiff(b[c[i,4],:], a[i,:])
k2 = setdiff(b[c[i,5],:], b[i,:])
c[i,6] = k1[1]
c[i,7] = k2[1]
end

return c
end

function test2(a,b)
c = hcat(a,zeros(Int64,size(a,1),4))

@inbounds @views for i = 1:size(a,1)
dex = Vector{Int64}()
for j = 1:size(b,1)
if a[i,:] ⊆ b[j,:]
dex = push!(dex,j)
end
end
c[i,4] = dex[1]
c[i,5] = dex[2]
end
# Determine the number of lines a belongs to b,
# here I am sure that there are only two lines in b containing a

@inbounds for i = 1:size(a,1)
k1 = setdiff(b[c[i,4],:], @view(a[i,:]))
k2 = setdiff(b[c[i,5],:], @view(b[i,:]))
c[i,6] = k1[1]
c[i,7] = k2[1]
end

return c
end


P.S.: I modified the code as suggested by @cchderrick, otherwise it will not work.

Then I get:

julia> a = [1 2 3; 4 5 6; 7 8 9];

julia> b = [1 2 3 5;
4 5 8 6;
7 0 9 8;
2 3 1 6;
4 5 6 9;
8 9 7 1];

julia> @btime test1($a,$b)
3.490 μs (105 allocations: 9.83 KiB)
3×7 Array{Int64,2}:
1  2  3  1  4  5  6
4  5  6  2  5  8  9
7  8  9  3  6  0  1

julia> @btime test2($a,$b)
2.318 μs (105 allocations: 7.20 KiB)
3×7 Array{Int64,2}:
1  2  3  1  4  5  6
4  5  6  2  5  8  9
7  8  9  3  6  0  1


and I am 100% sure that it can be improved. However, it will require knowledge of what the algorithm is doing.

2 Likes

yes!

If different iterations of your for loops are not related to each other use SIMD, to enhance the speed.
Check Matlab & Julia Matrix Operations Benchmark and particularly JuliaBenchSIMD for some inspiration on how to make use of SIMD. I have used SIMD in 3 MatrixAddition, MatrixMultiplication, ElementWiseOperations functions.

1 Like

Assuming that there are always two rows that match, no more and no less, and that size(a, 2) = 3 and size(b, 2) = 4, then I would replace the methods ⊆ and setdiff() by hard coded ones. You want to do this because if you do some “easy-profiling” like writing @time in front of the if a[i,:] ⊆ b[j,:] loops and in front of k1 = setdiff(b[c[i,4],:], a[i,:]) you will see that those are the lines where the script spends most of its time and almost all allocations. I also took the transpose the matrices a and b so that I’m accessing their elements column-wise (see the Performance Tips). @inbounds also helps.

My code is:

function foo2(a, b)
c = zeros(Int, 4, size(a, 2))
@inbounds for i=1:size(a, 2)
tmp = 1
for j=1:size(b, 2)
tmp > 2 && continue
is_1 = a[1, i] == b[1, j] || a[2, i] == b[1, j] || a[3, i] == b[1, j]
is_2 = a[1, i] == b[2, j] || a[2, i] == b[2, j] || a[3, i] == b[2, j]
is_3 = a[1, i] == b[3, j] || a[2, i] == b[3, j] || a[3, i] == b[3, j]
is_4 = a[1, i] == b[4, j] || a[2, i] == b[4, j] || a[3, i] == b[4, j]
# This line could potentially save time with larger matrix sizes.
#!is_1 & !is_2 & !is_3 & !is_4 && continue
is_123 = is_1 & is_2 & is_3
is_124 = is_1 & is_2 & is_4
is_134 = is_1 & is_3 & is_4
is_234 = is_2 & is_3 & is_4
if is_123
c[tmp, i] = j
c[tmp+2, i] = b[4, j]
tmp += 1
elseif is_124
c[tmp, i] = j
c[tmp+2, i] = b[3, j]
tmp += 1
elseif is_134
c[tmp, i] = j
c[tmp+2, i] = b[2, j]
tmp += 1
elseif is_234
c[tmp, i] = j
c[tmp+2, i] = b[1, j]
tmp += 1
end
end
end
return c
end

a = [1 2 3;
4 5 6;
7 8 9];

b = [1 2 3 5;
4 5 8 6;
7 0 9 8;
2 3 1 6;
4 5 6 9;
8 9 7 1];

aT = Matrix(transpose(a))
bT = Matrix(transpose(b))


Results:

>>> @btime foo($a,$b);
6.563 μs (103 allocations: 9.55 KiB)

>>> @btime foo2($a,$b);
173.424 ns (1 allocation: 176 bytes)


Nearly \times 40 faster. If you have problems with larger matrices then I’d try with multithreading in the outermost loop.

In any case, my script is not very elegant and is hard-coded for the input you gave, but maybe some of the above ideas help you.

1 Like

Maybe because I am not a native English speaker, there are some problems in expression, but your answer is what I want!!! Thank you very much .I hadn’t noticed before that “⊆” efficiency was very low.

### Here I would like to continue to ask:

• Are you talking about parallel processing at the outermost level like this: (because my actual matrix is really large)
using Distributed

function matrix_opt(a, b)
c = zeros(Int, 4, size(a, 2))
@inbounds @sync @distributed  for i=1:size(a, 2)
tmp = 1
for j=1:size(b, 2)
tmp > 2 && continue
is_1 = a[1, i] == b[1, j] || a[2, i] == b[1, j] || a[3, i] == b[1, j]
is_2 = a[1, i] == b[2, j] || a[2, i] == b[2, j] || a[3, i] == b[2, j]
is_3 = a[1, i] == b[3, j] || a[2, i] == b[3, j] || a[3, i] == b[3, j]
is_4 = a[1, i] == b[4, j] || a[2, i] == b[4, j] || a[3, i] == b[4, j]
# This line could potentially save time with larger matrix sizes.
#!is_1 & !is_2 & !is_3 & !is_4 && continue
is_123 = is_1 & is_2 & is_3
is_124 = is_1 & is_2 & is_4
is_134 = is_1 & is_3 & is_4
is_234 = is_2 & is_3 & is_4
if is_123
c[tmp, i] = j
c[tmp+2, i] = b[4, j]
tmp += 1
elseif is_124
c[tmp, i] = j
c[tmp+2, i] = b[3, j]
tmp += 1
elseif is_134
c[tmp, i] = j
c[tmp+2, i] = b[2, j]
tmp += 1
elseif is_234
c[tmp, i] = j
c[tmp+2, i] = b[1, j]
tmp += 1
end
end
end
return Matrix(transpose(c))
end

• I have another request similar to the title:
I have the following two matrices:
a=\left[ \begin{matrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ &...& \end{matrix} \right] , b=\left[ \begin{matrix} 1&4&3&2\\ 7&2&9&8\\ 4&6&1&5\\ 8&2&1&7\\ 11&3&4&8\\ &...&... \end{matrix} \right]

I am quite sure:

1. a are only three columns,b are only four columns.
2. In b there must be only one row containing a.For example:
the first line of a contains only the first line of b, and the second line of a contains only the third line of b.
3. I need to return the line number contained in a and b .And I also need to get a line in b that contains any elements other than a.So what I end up with is that c should be like this:
c=\left[ \begin{matrix} 1 & 2 & 3 &1&4\\ 4 & 5 & 6 &3&1\\ 7 & 8 & 9 &2&2\\ &&...& \end{matrix} \right]

I would like to ask you to help me write a piece of code in accordance with this requirement.
And I hope it execution efficiency is high (consider multi-core parallelism).

Thank you very much. I can see your idea, but my reading ability in the code is not very strong. I hope to get your help. Thank you.

Hi, I think it was not your English, but the code that was a little hard to get.

I would suggest @sync @distributed + SharedArrays or Threads.@threads in front of the for i=1:size(a, 2) loop. Please, try both and see which one works best for you.

I can see your idea, but my reading ability in the code is not very strong.

Many boolean expressions uncommented, sorry. The comments would be:

# True if the element b[2, j] is in a[:, i].
is_2 = a[1, i] == b[2, j] || a[2, i] == b[2, j] || a[3, i] == b[2, j]

...

# True if the elements of a[:, i] are contained in {b[1, j], b[3, j], b[4, j]}, no matter the order.
is_134 = is_1 & is_3 & is_4

...

# The element b[2, j] is not in a[:, i], so write b[2, j] in c.
elseif is_134
c[tmp, i] = j
c[tmp+2, i] = b[2, j]
tmp += 1


I would tell you to try to make the code you asked for now that I explained the code a bit more. The code you want is very similar to the previous one.

In case you still don’t know how to make it here it is

using BenchmarkTools

function foo(a, b)
c = zeros(Int, 2, size(a, 2))
@inbounds for i=1:size(a, 2)
for j=1:size(b, 2)
# True if the element b_{1, j} is in the ith row of a.
is_1 = a[1, i] == b[1, j] || a[2, i] == b[1, j] || a[3, i] == b[1, j]
# True if the element b_{2, j} is in the ith row of a.
is_2 = a[1, i] == b[2, j] || a[2, i] == b[2, j] || a[3, i] == b[2, j]
# True if ...
is_3 = a[1, i] == b[3, j] || a[2, i] == b[3, j] || a[3, i] == b[3, j]
is_4 = a[1, i] == b[4, j] || a[2, i] == b[4, j] || a[3, i] == b[4, j]
# This line could potentially save time with larger matrix sizes.
#!is_1 & !is_2 & !is_3 & !is_4 && continue
# True if the ith row of a is distributed in the first three elements of b[:, j]
is_123 = is_1 & is_2 & is_3
# True if the ith row of a is distributed in the elements: b[1, j], b[2, j], b[4, j].
is_124 = is_1 & is_2 & is_4
# True if the ith row of a is distributed in the elements: b[1, j], b[3, j], b[4, j].
is_134 = is_1 & is_3 & is_4
# True if ...
is_234 = is_2 & is_3 & is_4
if is_123
c[1, i] = j
c[2, i] = b[4, j]
elseif is_124
c[1, i] = j
c[2, i] = b[3, j]
elseif is_134
c[1, i] = j
c[2, i] = b[2, j]
elseif is_234
c[1, i] = j
c[2, i] = b[1, j]
end
end
end
return c
end

a = [1 2 3; 4 5 6; 7 8 9];
b = [1 4 3 2; 7 2 9 8; 4 6 1 5; 8 2 1 7; 11 3 4 8];
aT = Matrix(transpose(a))
bT = Matrix(transpose(b))

c = foo(aT, bT)
display(hcat(a, c'))

@btime foo($aT,$bT);
`

1 Like