0

can anyone help me with this?

Matlab Code

test = 300
f = @(u) 2+3*u-2*u.^3+u.^6-4*u.^8;
rng(12354)
u_test = unifrnd(-1,1,[test,1]);
y_test = f(u_test);
xt0 = 1*ones(test, 1);
xt1 = u_test;
xt2 = u_test.^2;
xt3 = u_test.^3;
xt4 = u_test.^4;
xt5 = u_test.^5;
xt6 = u_test.^6;
xt7 = u_test.^7;
xt8 = u_test.^8;
xt9 = u_test.^9;
x_test = [xt0 xt1 xt2 xt3 xt4 xt5 xt6 xt7 xt8 xt9];

theta_test = (x_test'*x_test)^-1*x_test'*y_test;

Python Code

import numpy as np
test = 3
min_test = -1
max_test = 1
f = lambda u: 2+3*u-2*u**3+u**6-4*u**8
np.random.seed(12345)
u_test = np.random.uniform(min_test, max_test, [test, 1])

y_test = f(u_test)

xt0 = np.ones([test, 1])
xt1 = u_test
xt2 = u_test ** 2
xt3 = u_test ** 3
xt4 = u_test ** 4
xt5 = u_test ** 5
xt6 = u_test ** 6
xt7 = u_test ** 7
xt8 = u_test ** 8
xt9 = u_test ** 9
x_test = np.array([xt0, xt1, xt2, xt3, xt4, xt5, xt6, xt7, xt8, xt9])

theta_test = (x_test.T * x_test) ** -1 * x_test.T * y_test

The final answer of theta in MATLAB code, which is also the correct answer, is a 10 * 1 matrix. It has 10 numbers in total, but the theta answer in the Python code is a 10 * 30 matrix. We have 300 numbers in the output.

How to do this multiplication with Numpy so that the final answer is the same as the MATLAB code answer?

mohamad halakoei
  • 3
  • 2
  • You've probably just missed a transpose in your Python implementation. That said, you're creating dynamic variable names, which is [bad, very bad](https://stackoverflow.com/a/32467170/5211833) for a plethora of reasons (both in MATLAB and Python). Simply use vector notation: `x_test = u_test.^(1:9)` – Adriaan Feb 01 '22 at 10:38
  • I'm having trouble running your MATLAB in Octave, but it looks like `x_test` is (n,10) size, and thus `(x_test'*x_test)^-1` should be (10,10). In MATLAB `*` is matrix multiplication. In `numpy` it is element-wise (`.*)` with `broadcasting`. `np.dot` or the matmul operator `@` is used for matrix multiplication. – hpaulj Feb 01 '22 at 17:20
  • There was a typo in MATLAB code that was fixed I also used np.dot, @ but to no avail – mohamad halakoei Feb 01 '22 at 17:33
  • My problem was with the initial stuff, generating the `u_test` and `y_test`. It would help a lot if you annotated both codes with the variables size/shape - intended and actual. – hpaulj Feb 01 '22 at 17:52

2 Answers2

0

As the comment on your question has said, it's better to write this code using an array rather than a list of variables.

That said, the issue is that you are using * which, in Python, denotes entrywise multiplication rather than matrix multiplication. What you should try instead is

theta_test = (x_test.T @ x_test) ** -1 * x_test.T @ y_test
Ben Grossmann
  • 4,387
  • 1
  • 12
  • 16
  • Unfortunately, I encountered this error....................... matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 3 is different from 10) – mohamad halakoei Feb 01 '22 at 17:30
  • Don't try to use `@` without first reading the `np.matmul` docs. And then carefully match the dimensions of your arrays with the function's specifications. – hpaulj Feb 01 '22 at 19:13
0

Your numpy code produces:

In [3]: x_test.shape
Out[3]: (10, 3, 1)
In [4]: y_test.shape
Out[4]: (3, 1)

Stripping off the last x_test dimension:

In [11]: x1 = x_test[:,:,0]
In [12]: (x1.T@x1).shape
Out[12]: (3, 3)
In [13]: np.linalg.inv(x1.T@x1)
Out[13]: 
array([[ 0.34307216, -0.63513713,  0.36347551],
       [-0.63513713,  8.44935183, -6.36062978],
       [ 0.36347551, -6.36062978,  5.43319377]])

Switching the transpose produces a (10,10) array, but it's singular

In [14]: np.linalg.inv(x1@x1.T)
Traceback (most recent call last):
  ...
LinAlgError: Singular matrix

The (10,3) can multiply the (3,1) to produce a (10,1)

In [18]: (x1@y_test).shape
Out[18]: (10, 1)

I'm having troubles running your MATLAB code in Octave, so can't generate its equivalent to see what you are trying to reproduce.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • There was a typo in MATLAB code that was fixed, you can take another look at it – mohamad halakoei Feb 01 '22 at 17:41
  • In translating code like this, I try to match the calculations step by step, especially getting the shapes to match. Don't just write a big script and hope the final results are right. Make a small test case, and track shapes and values step by step - in both languages. If you don't do that, I'll have to tag your question as missing debugging data. – hpaulj Feb 01 '22 at 18:00