9

I know how to solve A.X = B by least squares using Python:

Example:

A=[[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,0,0]]
B=[1,1,1,1,1]
X=numpy.linalg.lstsq(A, B)
print X[0]
# [  5.00000000e-01   5.00000000e-01  -1.66533454e-16  -1.11022302e-16]

But what about solving this same equation with a weight matrix not being Identity:

A.X = B (W)

Example:

A=[[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,0,0]]
B=[1,1,1,1,1]
W=[1,2,3,4,5]
double-beep
  • 5,031
  • 17
  • 33
  • 41
Eric H.
  • 2,152
  • 4
  • 22
  • 34
  • 1
    Have you looked at this link: http://stackoverflow.com/questions/19624997/understanding-scipys-least-square-function-with-irls – xnx Nov 25 '14 at 14:56
  • Yes; I tried: B=numpy.dot(B,W) before solving, but I have a message: numpy.linalg.linalg.LinAlgError: 0-dimensional array given. Array must be two-dimensional – Eric H. Nov 25 '14 at 15:03
  • 2
    If you take the dot product of two one-dimensional arrays you will get a scalar. Perhaps you mean to simply multiply the elements of B by those of W? Better to use numpy arrays rather than Python lists here. – xnx Nov 25 '14 at 15:08

3 Answers3

14

I found another approach (using W as a diagonal matrix, and matricial products) :

A=[[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,0,0]]
B = [1,1,1,1,1]
W = [1,2,3,4,5]
W = np.sqrt(np.diag(W))
Aw = np.dot(W,A)
Bw = np.dot(B,W)
X = np.linalg.lstsq(Aw, Bw)

Same values and same results.

Eric H.
  • 2,152
  • 4
  • 22
  • 34
  • 1
    the matrix product is more expensive than the element-wise product of @xnx's suggestion – DomTomCat Apr 03 '17 at 11:44
  • 1
    might be more computationally expensive but this is way more clear to read. +1 for code clarity with linear algebra – D A May 14 '17 at 23:46
  • Thanks! Why do you use sqrt(weight)? – MonsieurBeilto Oct 13 '18 at 23:22
  • 3
    @MonsieurBeilto, that's because in least squares method a sum of square displacements is minimized `(y - y0) ** 2`, therefore if you rescale y by `sqrt(w)`, a factor of `w` will pop up – icemtel Mar 29 '19 at 13:22
  • I am confused by the differences between this answer and xnx's answer. Why is the effect the same? You multiply the weights only to some entries of `A` while xnx multiplies them to all the entries of `A`, right? – lucidbrot Mar 05 '21 at 13:11
  • Nevermind, I got it. When multiplying the diagonal `W` with `A`, e.g. the entry at `[0,0]` is multiplied with every entry in the first row of `A` and all the other entries of the same column in `A` have no effect due to the zero-entries in `W`. It's simple matrix multiplication – lucidbrot Mar 05 '21 at 18:40
10

I don't know how you have defined your weights, but you could try this if appropriate:

import numpy as np
A=np.array([[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,1,1],[1,1,0,0]])
B = np.array([1,1,1,1,1])
W = np.array([1,2,3,4,5])
Aw = A * np.sqrt(W[:,np.newaxis])
Bw = B * np.sqrt(W)
X = np.linalg.lstsq(Aw, Bw)
xnx
  • 24,509
  • 11
  • 70
  • 109
  • 2
    Thanks. I don't know well numpy, and have to realize that W[:,np.newaxis] (or W[:,None]) gives a diagonal matrix. I.e. array([[1],[2],[3],[4],[5]]) means a 5x5 diagonal matrix with these values on the diagonal. – Eric H. Nov 26 '14 at 07:48
  • 1
    (Is it ? Because the result Aw is a 5x4 array, so that's the only explanation) – Eric H. Nov 26 '14 at 07:55
  • 2
    Sorry, I was confused by the * operator (I was not used to arrays). It is not a matrix product, but term to term. – Eric H. Nov 26 '14 at 09:08
  • 1
    That's right, NumPy arrays use * as the elementwise multiplication operator. You need to turn W into a column array to broadcast this multiplication correctly over A. There is also a NumPY matrix object, but it's more trouble than it's worth (IMO). – xnx Nov 26 '14 at 14:32
  • Thanks! Why do you use sqrt(weight)? Is it because it will be squared in error compute? – MonsieurBeilto Oct 13 '18 at 23:23
  • @MonsieurBeilto You can read about this in the Introduction section of the relevant Wikipedia page: https://en.wikipedia.org/wiki/Weighted_least_squares – if you assume the observational errors are uncorrelated then the weight matrix is diagonal with entries equal to the reciprocal of the measurement variances and the formulas simplify nicely if one refactors using the square root of this matrix. – xnx Oct 15 '18 at 07:05
0

scikit package offers weighted regression directly .. https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html#sklearn.linear_model.LinearRegression.fit

import numpy as np
# generate random data
N = 25
xp = [-5.0, 5.0]
x = np.random.uniform(xp[0],xp[1],(N,1))
e = 2*np.random.randn(N,1)
y = 2*x+e
w = np.ones(N)

# make the 3rd one outlier
y[2] += 30.0
w[2] = 0.0

from sklearn.linear_model import LinearRegression
# fit WLS using sample_weights
WLS = LinearRegression()
WLS.fit(x, y, sample_weight=w)

from matplotlib import pyplot as plt
plt.plot(x,y, '.')
plt.plot(xp, xp*WLS.coef_[0])
plt.show()

weighted regression without outlier

Karel Marik
  • 811
  • 8
  • 13