Векторизация создания матрицы последовательных степеней

Пусть x=1:100 и N=1:10. Я хотел бы создать матрицу x^N, чтобы столбец i th содержал записи [1 i i^2 ... i^N].

Я легко могу сделать это с помощью циклов for. Но есть ли способ сделать это с помощью векторизованного кода?


person alext87    schedule 13.10.2010    source источник
comment
+1 за возможность увидеть все эти разные подходы.   -  person zellus    schedule 14.10.2010


Ответы (5)


Я бы пошел на:

x = 1:100;
N = 1:10;
Solution = repmat(x,[length(N)+1 1]).^repmat(([0 N])',[1 length(x)]);

Другое решение (вероятно, гораздо более эффективное):

Solution = [ones(size(x)); cumprod(repmat(x,[length(N) 1]),1)];

Или даже:

 Solution = bsxfun(@power,x,[0 N]');

Надеюсь это поможет.

person Adrien    schedule 13.10.2010
comment
+1 для bsxfun + указатель на функцию + простые аргументы. Если вы можете выразить правило компактно, скорее всего, эта комбинация позволит вам выразить его в MATLAB ... - person Alex Feinman; 13.10.2010
comment
+1 за bsxfun () voodoo. @ Адриан, я думаю, это опечатка: N = 1:10; вместо N = 1: 0; - person zellus; 13.10.2010
comment
bsxfun () - самый приятный. Спасибо. - person alext87; 14.10.2010

Похоже на матрицу Вандермонда. Поэтому используйте vander:

A = vander(1:100);
A = A(1:10, :);
person Oliver Charlesworth    schedule 13.10.2010
comment
Единственная проблема заключается в том, что VANDER создает квадратные матрицы, поэтому вы выполняете довольно много дополнительной работы. - person gnovice; 13.10.2010
comment
@gnovice: не очень сложная проблема, чтобы ее обойти, см. код выше. - person Oliver Charlesworth; 13.10.2010
comment
Собственно, поскольку OP упомянул работу с квадратными матрицами в комментарии, эта функция выглядит хорошим выбором. Единственное изменение, которое вам нужно сделать, - это повернуть его с помощью ROT90, чтобы столбец i содержал степени i: A = rot90(vander(1:N)); - person gnovice; 14.10.2010

Поскольку ваши матрицы не такие большие, наиболее простой способ сделать это - использовать MESHGRID и поэлементный оператор мощности .^:

[x,N] = meshgrid(1:100,0:10);
x = x.^N;

Это создает матрицу 11 на 100, где каждый столбец i содержит [i^0; i^1; i^2; ... i^10].

person gnovice    schedule 13.10.2010
comment
хорошее использование сетки. Моя единственная проблема с этим методом и аналогичными, предложенными здесь, заключается в том, что они не используют i ^ -ю строку для вычисления строки i + 1 с небольшими затратами. Вариант cumprod делает это, но ему не хватает эстетики. - person Adrien; 14.10.2010
comment
@Adrien: Для больших матриц вы правы, что, вероятно, будет быстрее вычислять значения путем последовательного умножения на базовое значение. Для относительно небольших матриц, подобных приведенному здесь примеру, кажется разумным компромиссом пойти на что-то легкое для чтения и понимания, но немного менее эффективное. ;) - person gnovice; 14.10.2010
comment
Мне нравится этот метод ... В моем реальном программировании я создаю 10000 на 10000. Так что, возможно, мне следует использовать более эффективные подходы. Тем не менее, это полезно знать. Спасибо! - person alext87; 14.10.2010
comment
@ alext87: Это очень большие матрицы, и вы, вероятно, столкнетесь с проблемы с памятью. Фактически, они настолько велики, что когда я создаю матрицу такого размера, памяти не остается даже для выполнения с ней базовых операций. Например, clear all; V=ones(10000); V=V.'; выдает ошибку нехватки памяти. Кроме того, нет никакого смысла увеличивать число до такой большой степени. Даже просто 2^1000 приводит к 1.0715e+301, что немного меньше максимального значения для удвоений. - person gnovice; 14.10.2010
comment
@Adrien: На самом деле функция VANDER, которая была предложенный Оли Чарльзуортом использует преимущества упомянутых вами эффективных вычислений. Вы можете увидеть код, набрав type vander в командном окне MATLAB. - person gnovice; 14.10.2010

Не уверен, что это действительно соответствует вашему вопросу.

bsxfun(@power, cumsum(ones(100,10),2), cumsum(ones(100,10),1))

РЕДАКТИРОВАТЬ: Как указал Адриан, моя первая попытка не соответствовала вопросу OP.

xn = 100;
N=10;
solution = [ones(1,xn); bsxfun(@power, cumsum(ones(N,xn),2), cumsum(ones(N,xn),1))];
person zellus    schedule 13.10.2010
comment
Первая строка матрицы вычислений должна быть одна, в вашем случае это 1:100 - person Adrien; 13.10.2010

Почему бы не использовать простой для понимания цикл for?

c = [1:10]'; %count to 100 for full scale problem
for i = 1:4; %loop to 10 for full scale problem
    M(:,i) = c.^(i-1)
end

Чтобы понять умные векторизованные версии этого кода, которые показали люди, нужно больше думать. Мой способ ведения дел скорее варварский, но любой, кто его прочитает, поймет это.

Я предпочитаю простой для понимания код.

(да, я мог бы выделить заранее. Не стоит понижать ясность для таких небольших случаев.)

person MatlabDoug    schedule 13.10.2010
comment
Если вы используете цикл, вы должны заполнять матрицу построчно, а не использовать оператор. ^. Вместо этого заполните первую строку единицами, а затем создавайте каждую новую строку, умножая предыдущую на [1:100]. Тогда у вас будет и удобочитаемость, и эффективность. - person Adrien; 14.10.2010
comment
Предложение Адриана увеличивает скорость на 12%. Для обсуждаемого здесь масштаба проблемы это было незаметно. Мне пришлось запустить это 100000 раз, чтобы все это заняло 2,8 против 2,5 секунды. На мой взгляд, мою версию проще читать, так как для M нет шага инициации, но и не о вкусе! :) - person MatlabDoug; 15.10.2010