Вы можете создать одну систему уравнений, содержащую множество независимых «подсистем» уравнений, поместив подматрицы a1 по диагонали матрицы 12000x12000 следующим образом:

a1(1,1) a1(1,2) a1(1,3)    0       0      0        0       0      0       
a1(2,1) a1(2,2) a1(2,3)    0       0      0        0       0      0       
a1(3,1) a1(3,2) a1(3,3)    0       0      0        0       0      0       
   0       0       0    a1(1,4) a1(1,5) a1(1,6)    0       0      0       
   0       0       0    a1(2,4) a1(2,5) a1(2,6)    0       0      0       
   0       0       0    a1(3,4) a1(3,5) a1(3,6)    0       0      0       
   0       0       0       0       0       0    a1(1,7) a1(1,8) a1(1,9)   
   0       0       0       0       0       0    a1(2,7) a1(2,8) a1(2,9)   
   0       0       0       0       0       0    a1(3,7) a1(3,8) a1(3,9)

а затем разделите его на a2(:).

Это можно сделать с помощью kron и разреженной матрицы, подобной этой (источник):

a1_kron = kron(speye(12000/3),ones(3));
a1_kron(logical(a1_kron)) = a1(:);
a = a1_kron\a2(:);
a = reshape(a, [3 12000/3]);

Преимущество — скорость: это примерно в 3-4 раза быстрее, чем цикл for с предварительным выделением на моем ПК.

Недостаток. При таком подходе следует учитывать один недостаток: при использовании левого деления Matlab ищет наилучший способ решения систем линейных уравнений, поэтому, если вы решаете каждую подсистему независимо, для каждой подсистемы будет выбран лучший путь, но если вы решаете тему как одну систему, Matlab найдет лучший путь для всех подсистем вместе, а не лучший для каждой подсистемы.

Примечание. Как показано в Stefano M ответ, использование одной большой системы уравнений (с использованием kron и разреженной матрицы) быстрее, чем использование цикла for (с предварительным выделением) только для очень небольшого размера под- системы уравнений (на моем компьютере число уравнений ‹= 7) для больших размеров подсистем уравнений использование цикла for быстрее.

Сравнение различных методов

Я написал и запустил код для сравнения 4 разных методов решения этой проблемы:

1. для цикла, без предварительного распределения
2. для цикла, с предварительным выделением
3. kron
4. cellfun

Контрольная работа:

n = 1200000;
a1 = rand(3,n);
a2 = rand(3,n/3);

disp('Method 1: for loop, no preallocation')
tic
a_method1 = [];
for ii = 1:3:n
  a_method1 = [a_method1 a1(:,ii:ii+2)\a2(:, ceil(ii/3))];
end
toc
disp(' ')

disp('Method 2: for loop, with preallocation')
tic
a1_reshape = reshape(a1, 3, 3, []);
a_method2 = zeros(size(a2));
for i = 1:size(a1_reshape,3)
    a_method2(:,i) = a1_reshape(:,:,i) \ a2(:,i);
end
toc
disp(' ')

disp('Method 3: kron')
tic
a1_kron = kron(speye(n/3),ones(3));
a1_kron(logical(a1_kron)) = a1(:);
a_method3 = a1_kron\a2(:);
a_method3 = reshape(a_method3, [3 n/3]);
toc
disp(' ')

disp('Method 4: cellfun')
tic
a1_cells = mat2cell(a1, size(a1, 1), repmat(3 ,1,size(a1, 2)/3));
a2_cells = mat2cell(a2, size(a2, 1), ones(1,size(a2, 2)));
a_cells = cellfun(@(x, y) x\y, a1_cells, a2_cells, 'UniformOutput', 0);
a_method4 = cell2mat(a_cells);
toc
disp(' ')

Полученные результаты:

Method 1: for loop, no preallocation
Elapsed time is 747.635280 seconds.
 
Method 2: for loop, with preallocation
Elapsed time is 1.426560 seconds.
 
Method 3: kron
Elapsed time is 0.357458 seconds.
 
Method 4: cellfun
Elapsed time is 3.390576 seconds.

Сравнив результаты четырех методов, можно увидеть, что использование метода 3 — kron дает немного разные результаты:

disp(['sumabs(a_method1(:) - a_method2(:)): ' num2str(sumabs(a_method1(:)-a_method2(:)))])
disp(['sumabs(a_method1(:) - a_method3(:)): ' num2str(sumabs(a_method1(:)-a_method3(:)))])
disp(['sumabs(a_method1(:) - a_method4(:)): ' num2str(sumabs(a_method1(:)-a_method4(:)))])

Результат:

sumabs(a_method1(:) - a_method2(:)): 0
sumabs(a_method1(:) - a_method3(:)): 8.9793e-05
sumabs(a_method1(:) - a_method4(:)): 0

Вы решаете серию из N систем с m линейными уравнениями в каждой, N систем имеют вид

Ax = b

Вы можете преобразовать их в единую систему линейных уравнений Nm:

|A1 0  0  ... 0 | |x1|   |b1|
|0  A2 0  ... 0 | |x2|   |b2|
|0  0  A3 ... 0 | |x3| = |b3|
|.  .  .  ... . | |. |   |. |
|0  0  0  ... AN| |xN|   |bN|

Однако решение этой одной системы уравнений намного дороже, чем решение всех маленьких. Как правило, стоимость составляет O(n^3), поэтому вы переходите от O(Nm^3) к O((Nm)^3). Огромная пессимизация. (Элиаху доказал, что я ошибаюсь здесь, по-видимому, можно использовать разреженность матрицы.)

Сократить вычислительные затраты можно, но необходимо предоставить гарантии относительно данных. Например, если матрицы A положительно определенные, системы могут быть решены дешевле. Тем не менее, учитывая, что вы имеете дело с матрицами 3x3, выигрыш в них будет невелик, так как это довольно простые для решения системы.

Если вы спрашиваете об этом, потому что считаете, что циклы в MATLAB по своей природе медленные, вы должны знать, что это уже не так, и не было так с тех пор, как MATLAB получил JIT 15 лет назад или около того. Сегодня многие попытки векторизации приводят к столь же быстрому коду, а часто и к более медленному коду, особенно для больших данных. (Я мог бы найти некоторые тайминги, которые я разместил здесь, на SO, чтобы доказать это, если это необходимо.)

Я думаю, что решение всех систем за один раз может уменьшить количество проверок, которые MATLAB выполняет каждый раз, когда вызывается оператор \. То есть жесткое кодирование размера и типа проблемы может улучшиться во всем. Но единственный способ сделать это — написать MEX-файл.

ПредельноСамым большим улучшением будет предварительное выделение выходной матрицы вместо ее увеличения:

A1 = reshape(A1, 3, 3, []);
a = zeros(size(A2));
for i = 1:size(A1,3) 
    a(:,i) = A1(:,:,i) \ A2(:,i);
end

С предварительно распределенным массивом, если доступен Parallel Toolbox, вы можете попробовать parfor

Решение разреженной матрицы другого ответа Элиаху Аарона действительно очень умно, но его преимущество в скорости не является общим, а зависит от конкретного размер проблемы.

С помощью этой функции вы можете исследовать различные размеры задач:

function [t2, t3] = sotest(n, n2)

a1 = rand(n,n*n2);
a2 = rand(n,n2);

tic
a1_reshape = reshape(a1, n, n, []);
a_method2 = zeros(size(a2));
for i = 1:size(a1_reshape,3)
    a_method2(:,i) = a1_reshape(:,:,i) \ a2(:,i);
end
t2 = toc;

tic
a1_kron = kron(speye(n2),ones(n));
a1_kron(logical(a1_kron)) = a1(:);
a_method3 = a1_kron\a2(:);
a_method3 = reshape(a_method3, [n n2]);
t3 = toc;

assert ( norm(a_method2 - a_method3, 1) / norm(a_method2, 1) < 1e-8) 

Действительно, для n=3 метод разреженных матриц явно лучше, но для увеличения n он становится менее конкурентоспособным.

Приведенный выше рисунок был получен с

>> for i=1:20; [t(i,1), t(i,2)] = sotest(i, 50000); end
>> loglog(1:20, t, '*-')

Мой последний комментарий заключается в том, что явный цикл с плотным оператором \ действительно быстр; формулировка разреженной матрицы немного менее точна и может стать проблематичной в крайних случаях; и наверняка решение разреженной матрицы не очень читабельно. Если количество n2 систем, которые нужно решить, очень-очень велико (>1e6), возможно, следует изучить специальные решения.